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SECTION  I 
INTRODUCTION 


This  report  is  concerned  with  theories,  formulations,  and  computer 
programs  for  the  analysis  of  shell  structures  in  general  and  fiber  rein¬ 
forced  rocket  motor  cases  in  particular.  Such  analysis  is  intended  for 
use  in  efficient  design  of  missile  structures  consisting  of  complex 
structural  components  subjected  to  various  loading  conditions.  Only  the 
isothermal  problems  ore  reported  herein. 

Various  material  properties  such  as  elasticity,  plasticity,  visco¬ 
elasticity,  and  viscoelastoplasticity  are  considered  in  the  present  study. 
Schapery  Cl9,13],  among  others,  confirmed  that  glass  fiber-reinfor¬ 
ced  structures  embedded  in  epoxy  exhibit  significant  rate-dependent  elasto 
plastic  or  viscoelastoplastic  behavior. 

The  response  of  such  structures  depends  greatly  on  types  of  loading, 
static  or  dynamic.  Perzyna  [18]  presented  the  thermodynamic  foundations 
of  the  theory  of  viscoplasticity.  Experimental  evidet.ces  in  his  study 
necessitated  the  simultaneous  consideration  of  viscoelastic  and  plastic 
properties  of  e  material.  Both  static  and  dynamic  loadings  will  be  stu¬ 
died  in  detail  in  the  present  study. 

Yielding  of  the  fiber-reinforced  structure  in  the  context  of  present 
state  of  art  still  remains  a  controversial  matter.  Studies  of  plastically 
anisotropic  materials  were  initiated  by  Hill  [8  ]  who  postulated  the  form 
of  yield  condition  based  on  the  von  Mlses  criterion  for  isotropic  plastic 
material.  Hu  [lO],  on  the  other  hand,  generalized  the  Tresca  shear-stress 
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yield  criterion  in  a  similar  manner.  These  theories,  however,  are  incap¬ 
able  of  taking  into  account  the  presence  of  fibers  in  a  ductile  matrix. 
Mulhern,  Rogers  and  Spencer  [  14  ]  proposed  a  general  continuum  theory  of 
the  mechanical  behavior  of  fibers  which  are  considered  plastically  inex- 
tensible.  Lance  and  Robinson  Q.2  ]  more  recently  presented  a  theory  of 
ductile  behavior  of  composite  materials  based  on  physical  ideas  related  to 

those  contained  in  the  maximum  shear  stress  theory  of  isotropic  materials. 

/ 

It  is  evident  that  this  later  theory  removes  shortcomings  of  all  other 
theories  but  successful  application  to  a  shell  structure  with  fiber  angle 
plys  appears  to  be  doubtful.  In  view  of  these  developments,  the  most  pro¬ 
mising  approach  is  to  modify  Hill’s  theory  to  incorporate  strain-hardening 
of  fibers  in  similar  context  of  the  work  of  Lance  and  Robinson.  Jensen 
[11]  and  Whang  [23]  applied  this  idea  to  orthotropic  material.  In  the 
present  study  a  further  extension  is  made  to  account  for  layers  of  fiber 
angle  plys. 

Thin  and  thick  shells  are  considered.  A  treatment 
as  a  thin  shell  requires  plane  stress  approximation  of  the  finite  element 
model  in  which  anisotropic  parameters  for  each  axisymmetric  layer  of 
angle  plys  through  the  shell  thickness  must  be  characterized.  Detailed 
study  on  this  subject  which  has  not  appeared  in  the  literature  will  be 
elaborated  in  the  present  report.  The  results  of  a  thin  shell  approxi¬ 
mation  are  compared  with  three  dimensional  thick  shell  theory  in  which 
imposition  of  axisymmetry  permits  use  of  a  plane  strain  isoparametric 
finite  element.  This  comparison  points  to  a  distinct  superiority  of  three 
dimensional  theory  in  dealing  with  filament  wound  fibers.  The  computer 
program  is  capable  of  analyzing  a  body  composed  of  a  combination  of  several 
monolithic  materials  and  composice  materials. 
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Based  on  these  ideas  and  the  theoretical  background,  details  of  form- 
filiation  of  governing  aquations  are  presented  in  the  following  sections.  A 
Chin  shell  is  discussed  in  Section  2  and  a  thick  shell  in  Section  3  with 
specialized  topics  included  in  Appendices  A.l  through  A. 12.  Computer 
programs  are  written  with  maximum  utilization  of  currently  available  numer¬ 
ical  methods.  Demonstrative  problems  are  solved  and  results  presented  in 
Section  4.  Conclusions  and  recommendations  are  given  in  Section  5.  Docu¬ 
mentation  of  the  computer  programs  is  included  in  Appendices  B.l  through 
B.6. 
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SECTION  2 


THIN  SHELL  ANALYSIS 


2.1  GENERAL 

Since  an  extensive  review  of  developments  of  thin  shell  theory 
is  beyond  the  scope  of  the  present  study  we  simply  begin  with 
Novozhilov's  theory  [15,16]  which  is  relatively  widely  accepted  among 
the  engineers.  Modifications  a  re  then  Introduced  to  account  for 
anisotropy,  fiber-reinforced  composites,  plasticity,  viscoelasticity, 
and  viscoelastoplasticity.  Axisymmetric  shells  under  static  and  dy¬ 
namic  loadings  will  be  considered. 

2.2  CONSTITUTIVE  EQUATIONS  FOR  ELASTIC  AND  PLASTIC  BEHAVIOR 
Based  on  the  Love-Kirchhof f  hypothesis,  small  strains,  large 

displacements,  and  moderate  rotations,  the  membrane  strain  tensor 
e  _  and  the  bending  strain  tensor  v  are  given  below: 


2  fu„|P  +  "ri,,  -  +  <u:<v  +  "Vt*:**  “Vvi 


(2.1) 


I  ft  "  Un|«  b*  “  u.bX|a  "  u  I  +  b«  b  u3 

Wp  X  B  o/  )  nr  A  P  u  Irr  f  unt 


(2.2) 


Here  all  Greek  letters  range  from  1  to  2  and  u  and  u'  are  the  dis¬ 
ci t 

placements,  b  is  the  second  fundamental  tensor  ,  the  commas  and 
strokes  represent,  respectively,  ordinary  and  covariant  differentia¬ 


's. 


tions.  It  should  be  noted  that  all  nonlinear  terms  in  the  bending 


strain  are  neglected.  Derivations  of  (2.1)  and  (2.2)  and  the  phys¬ 
ical  components  based  on  Novozhilov's  theory  are  given  in  Appen¬ 
dix  A.I. 

For  linear  elastic  behavior,  the  stress  tensor  is  given  by 


O' A  _  'Yp  .  0<P 

<1  "<,}  +  f>(b)  <2*3> 

where  the  membrane  stress  tensor  off)  and  the  bending  stress  tensor 
a?J)  are,  respectively, 


CYR 

<*<■> 


h#PXue 

A|l 


<y|l 

<*<b> 


hfl  n  >PXll 
12  D 


X('. 


(2.4) 


in  which  h  is  the  shell  thickness  and  is  the  tensor  of  elastic 

moduli. 

If  yielding  of  the  material  is  considered  we  must  modify  (2.4) 
in  order  to  incorporate  plastic  behavior.  For  the  isotropic  von 
Mises  material  the  plastic  potential  function  f  becomes  [8J 


f  -  -  3J  (2.5) 

where  J  is  the  second  deviatoric  stress  invariant  and  a  is  the 
equivalent  yield  stress.  For  the  anisotropic  von  Mises  material 
[8l>  however,  we  have 

£  ’  »  *  I  *,  ,»/  S’  ’s"  (2-6) 


where  is  th«i  anisotropic  parameters,  S* J  is  the  deviatoric 
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stress  tensor,  and  *-  1,2,3.  For  the  case  of  plane  stress 

(2.6)  becomes 


f  ■  ?  ’  |  aTA  a 

in  which 


and 


O  -  [  a11  aa?  | 


0  0  «V,J 


(2.7) 


Explicit  forms  of  (2.5)  and  (2.7)  are  given  in  Appendix  A.2. 

The  incremental  plastic  strain  tensor  is  related  by  an  asso¬ 
ciated  flow  rule. 


where  d\  is  the  positive  constant. 

Differentiating  (2.5)  or  (2.7)  yields 


> 


where 


2 ^ 


(2.8) 


(2.9) 


Let  the  plastic  modulus  be  given  by 


7 


dy 


(2.10) 


where  dy  Is  Che  incremental  equivalent  yield  strain.  By  equating 

(  p) 

the  incremental  plastic  work  dW  done  by  the  current  stress  and 
current  incremental  plastic  strain  with  that  done  by  the  equivalent 
yidld  stress  and  incremental  equivalent  yield  strain, 


(p)  B  f  p  )  p  ) 

dW  -  cr  dy  B  *  crdY 
Of  “ 


(2.11) 


and  substituting  (2.8)  into  (2.11),  we  obtain 


*•(  p) 

d\  ■  dY  /2tr 


(2.12) 


Now  Inserting  (2.12)  into  (2.8)  yields 


<p>  -<p> 

dY  I  ■  Z 

cv*  nrp  arp 


(2.13) 


Since  the  total  incremental  strain  is  the  sum  of  the  elastic 
and  plastic  strains  we  can  write 


(.)  (p) 

dy  .  -  dy  .  -  dY  . 

T<*P  aP  orP 


(2.14) 


where  dy  .  is  the  elastic  strain  tensor.  This  enables  us  to  ex- 

«l» 

press  the  total  Incremental  stress  in  the  form, 


•  -A;’) 


(2.15) 


2.3  CONSTITUTIVE  EQUATIONS  FOR  VISCOELASTOPLASTIC  BEHAVIOR 

The  principle  of  conservation  of  energy  states  that  the  time 
rates  of  the  kinetic  energy  k  and  the  internal  energy  U  are  equal 
to  the  mechanical  power  R  and  the  heat  energy  Q.  In  the  absence 
of  the  thermal  loading  this  principle  may  be  expressed  as 


J 


it  +  U  -  R 


(2.20) 
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Here  the  superposed  dot  represents  the  time  rate  and 


«  1  f  . 

-  r  I  pv  v  4 
*  I  nt  at 


»  If 


(2.21) 


(2.22) 


R-  pPvdV  +  l  ."Nv 


(2.23a) 


in  which  p  is  the  density;  v  is  the  velocity  component;  $  is  the 

a 

internal  energy  density;  is  the  body  force;  is  the  surface 
traction;  and  n  is  the  unit  normal  to  the  surface.  Here  the  small 

(X 

strain  and  rectangular  cartesian  coordinates  are  used.  Using  the 
Green*- Gauss  theorem,  (2.23a)  becomes 


R  *  I  (pF\  +  v  +  a*  v  )dV 


(2.23b) 


Now,  Inserting  (2.21)  and  (2.23b)  into  (2.20)  yields 


f  £  a<xJ  +  ’  na*>Vp  ’  +  aaPvR,Jdv "  ° 


(2.24) 


For  the  principle  of  linear  momentum  to  hold  and  for  arbitrary 
volumes  we  must  have 


aa*  +  nF*  -  pa"  -  0 
>o- 
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> 


and 

> 

a 


(2.25) 


Here  the  commas  denote  ordinary  differentiation  and  aa  is  the 


acceleration. 

Our  objective  here  is  to  propose  a  form  of  free  energy  func¬ 
tions  in  incremental  quantity  such  that  the  non-smooth  or  inelastic 
strains  may  be  included  for  a  small  time  interval  At.  For  iso¬ 


thermal  conditions,  the  incremental  free  energy  <s(At)  and  stresses 


ad 

a 


(At)  are  assumed  to  be  functions  of  incremental  strains  V  _ 

(V  p 


(At)  - 


+  v,U>(At> 

(  r) 

variables)  o',  j  (At) 


and  incremental  internal  variables  (hidden 


(rX.V  . 

*  or,  j  (At)  + 


( r  X  p ) 

j  (At). 


This  statement  may 


be  written  as 


( « ) 


<pi 


♦  (At)  •♦[Y,J(At),  Y,  j  (At)  ,  ,vn  (At), 


(  r  )  (  p  ) 
‘  J 


<*« .  (At) | 


(2.26) 


<p>. 


o‘i(At)  -  S[Yu(At),  Yt  j  (At), 


(r)(i)  .  (r)<Pl 


<vx  j  (At), 


Of 


i  i 


(At)] 


(2.27) 


For  isothermal  conditions,  the  free  energy  is  the  same  as  the 
internal  energy  so  that 

o*  -  pi  ■ 

or  for  the  small  time  interval  At, 
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pi  (At)  -  oa9(bt)  {Y^IaO  +  V^(At)-} 


(2.28) 


At  this  point  we  introduce  hero  the  incremental  form  of  free 
energy  in  a  truncated  Taylor  series  expansion  [5  ], 


C*8  Xu 


V  )(%  +  V  1 


n 

V  <r>(P)w  <•>  <p\ 

z  o*<  V  +  v  K\»  + V,1 


(2.29) 


wliere  are  stiffness  constants  associated  with  the  internal 

variables.  Note  that  (2.29)  has  the  form  of  truncated  Taylor  series 

expansion  only  to  include  quadratic  terms.  However,  the  product 

term  of  Y* *  *  and  Y  **  is  missing.  This  is  because  an  explicit  mater- 

i*P  <*P 

(•)  (p) 

ial  kernel  relating  the  product  of  y  and  Y  is  nonexistent1  and 

a* 

coupling  of  elastic  and  plastic  strains  can  be  obtained  using  any 

(  T  ) 

one  of  the  failure  theories.  Lastly,  q  defined  here  as  the  iriter- 

QfP 

nal  variables  represent  time  dependent  physico-chemical  properties 
or  simply  a  viscous  behavior  which  may  be  expressed  as 


V  */ 


■<C'T)  .  ,  .  . 

“5 -  Y  .  (  r)d  t 


(2.30) 


where  r  is  the  time  variable  and  T(  j  is  the  relaxation  time.  In 
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order  to  facilitate  an  explicit  integration  we  assume  a  linear 

variation  of  within  the  time  interval  At  given  by 
arp 


*>)-  -  y8)  > 


(2.31) 


where  (si  is  the  current  time  step.  Substituting  (2.31)  in  (2.30) 
and  performing  integration  we  obtain 


(f)  <n(r)  ,  v  (r).  ,  .  (f). 

%,f>  ■  *  V"l>  +  8  V„,M  +  C  vc,1.{s> 


(2.32) 


in  which 


<t>  ,  -At 

A  *  exp(-= -  )  , 


(f)  (r>  (  r  ) 

B  -  T(r)(  i  -  A  ) 


(  r  )  (  r  ) 

P  «  T(r)(l-  *  ) 


(r)  (r) 

t  ■  TT~(1"  A  > 


The  derivation  of  these  parameters  is  given  in  Appendix  a«3. 
Rewriting  (2.28)  for  the  current  time  step  (s)  as 


f ***&+%&  XV i 


+  V  9\, 


a  p(*)  (Y  W  +  Yap(«0)  =0 


(2.33) 


and  substituting  (2,32)  and  (2.29)  into  (2.33)  yields 


r-l 


(  r  )  (  «  } 

B  V  »- 

X(.«. 


D  + 


(r>.  , 

C  V,  fe) ) 

All 


(  T  )  • 

+  q  <S)  V  CS)  + 
Xii  arp 


<S)  . 

v  flq,  cs>) 

ap  Xu 


v‘ ’  ®  *6*nltri?o>  £?<» 


.(•) 

Since  all  variations  other  than  Y  are  not  arbitrary  we  must  have 

crp 

the  relationship 


**18)  -  ('a'^’cs-D  +C6^  cs-1)  +<C>  v  *  <s> ) 

X(i  *(r)  Xti  X|i  Xp, 


Xpt 

(2.34) 


and 


( f  )  (  *  )  (  p  )  (p)(p)#(p) 

-  -  %<*> v» 


(r) 


+  q  Cs)  ^  ca)  +  V  CO)  $  (s)  j  *  0 
Qtf  Xu  Qffl  \u 


(2.35) 


It  should  be  noted  that  (2.34)  results  from 


a  p  T  6* 
a  p  *  P— ) 

BY  . 

ap 


which  states  that  the  stresses  are  derivable  from  the  free  energy 


functions.  In  our  specific  problem,  however,  this  stress  is  due  to 
an  elastic  strain  and  a  law  governing  the  plastic  strain  is  needed 
to  obtain  the  stress  due  to  a  total  strain  as  already  demonstrated 
in  Section  2.2.  The  relationship  (2.35)  may  be  considered  as  the 
dissipation  which  plays  a  significant  role  in  heat  conduction  prob¬ 
lems.  However,  for  the  isothermal  problems  as  considered  here  these 
terms  of  (2.35)  are  not  needed  in  the  analysis. 

Now,  if  the  plastic  behavior  is  considered  we  must  calculate 
the  stress  incrementally.  For  static  analysis  the  load  is  applied 
in  increments.  Far  transient  analysis  or  time  dependent  problems, 
the  total  loada  are  applied  for  each  time  increment.  In  either  case 
the  incremental  stresses  da°^  are  calculated  iteratively  until  con¬ 
vergence  is  achieved.  Therefore,  we  must  express  (2.34)  in  incre¬ 
mental  form  associated  with  the  total  strain  in  order  to  be  combined 
with  (2.18)  for  iterative  cycling  upon  direct  superposition  of  visco¬ 
elastic  behavior  and  plastic  behavior.  Such  superposition  is  achieved 
*:y  calculating  (jP**11'  initially  by  viscoelastic  stresses.  These  ar¬ 
guments  require  that  the  second  term  of  (2.18)  for  the  current  time 
step  (s)  la  added  to  the  incremental  form  of  (2.34), 

do0^  -  +V  ^Adqfs-l)  +(B)dV(s-l> 

Xu  ( r )  X(i  X(i 

+  C  iy<B)  ]  +  fP^dYf©  (2.36) 

Xu  Xu 

In  earlier  studies,  the  authors  obtained  a  slightly  different  form 


M. 


for  drj*%>  by  means  of  the  generalized  Maxwell  model  [1].  However, 
a  proper  choice  of  the  relaxation  time  would  yield  identical  results. 


2.4  THE  FINITE  ELEMENT  EQUATIONS  OF  MOTION 

The  use  of  the  finite  element  method  is  widespread  [17,25].  The 
advantage  is  due  to  its  capability  in  handling  complex  geometries, 
boundary  and  loading  conditions. 

The  generalized  coordinates  at  o  node  placed  in  the  midsur¬ 
face  of  the  ring  element  include  the  displacements  in  the  meridi¬ 
onal,  tangential  and  transverse  directions  plus  the  meridional 
directions  .The  tangential  displacement  need  not  be  considered  for 
axisymmetrical  deformations.  For  future  use  in  connection  with 
asymmetrical  loadings  and  for  the  purpose  of  generality,  we  shall 
Include  the  tangential  displacement  in  our  formulation.  The  element 
generalized  coordinates  are  given  by 


0,  -  *1N0N 


(2.37) 


where  ^tM  is  the  normalized  displacement  functions  and  the  super¬ 
script  ranges  from  1  to  the  total  number  of  generalized  coordinates. 
For  the  two  node  element  as  in  the  case  of  the  meridional  line  ele¬ 
ment  used  in  the  present  study  we  have  n*  2t  for  the  maximum  range. 

Appendix'  A.4  gives  the  explicit  form  of  ♦  for  linear  variations 
of  meridional  and  tangential  displacements,  cubic  variation  of  trans¬ 
verse  displacement  and  the  corresponding  meridional  rotation.  Such  an 
element  has  proven  to  be  quite  efficient  [24]. 
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j 


The  general  finite  element  equation  of  motion  is  obtained  by 
using  (2.3),  (2.37),  (2.25),  (2.22)  and  (2.33b)  in  (2.20), 


(2.38a) 


(2.38b) 


Here  Hnh  is  the  mass  matrix 

p*‘*1HdV  (2.39) 

and  PN  is  the  nodal  generalized  force  derived  from  (2.23b).  It 
should  be  noted  that  (2.38b)  is  obtained  from  the  relationship 


Y  ,  ■  e  +S  X  _  ,n  /  n  „\ 

c*P  afl  (2.40a) 

and  the  Integration  through  the  thickness  coordinate  £  is  contained 
in  ogj  and  of  (2.38b).  In  view  of  (2.1),  (2.2),  and  (2.37), 

we  have 


e  .  *  Am  e  +  a©*«5N 


(2.40b) 
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and 


.  m 

M«  9  N*P  +  CN«(y||0  '  f0H  “  ®NaP 


where  AMflfp,  CM(ta|j  ,  and  BNCyp  represent  derivatives  of  the  normalized 
interpolation  function  in  the  membrane  and  bending  strains  of  (2.1) 
and  (2.2). 

We  now  introduce  a  perturbation  or  a  variation  to  (2.38b)  for 
an  arbitrary  time  step  (s) 


Msad«%>  + 


'f*'a<(!)B>  <A"aP 


+  c, 


nm» 


)dA  +  jc^c 


Cs)  CNHapdA  dQ(B) 


pdA  «  dPN(S) 


where  d<^*^>  and  dO^P^ooare  deduced  from  (2.36), 


(2.41) 


n 

<n#Bv  -  haf**  +  fr^)de^+ h  ^;vic<£)dys) 


+  d^q^'s-l)  +  ^di^tP-1)] 


(2.42a) 


I 


and 


*o?!p> 


-  ^Xu  + 


dxc 


(C>dx<9) 


,  ({)>)(b)  (r)  . 

+  Ad^  (s-l)  +  6  dx«-l)  (2.42b) 


1 
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Substituting  (2.42a)  and  (2.42b)  Into  (2.41)  the  local  finite  element 
equation  of  motion  becomes 


+  JM„d<te)+  Kn*  dQ(jj)  - 


(2.43) 


in  which  JNN  and  KNM  are  the  viscosity  matrix  and  linear  elastic 


stiffness  matrix,  respectively, 


JN*  "  h 


VsW* +  nf  2 

tm\  JA  rml 


bn<*Pbh}^a 


|  +  12  dA 

A  A 


(2.44) 


(2.45) 


The  nodal  vector  dF^sj  contains  not  only  the  external  load  but  also 
the  pseudo  loads  of  various  sources , 


,  (V)  (8)  (p)  .  (N) 

dF^EF)  ■  dPN(p)  +  dFN  tsj  +  dFN  is)  +  dFN  is)  +  dFN  (s) 


(V)  (8)  h  ( p )  M  (  H ) 

dF^)  -  dP^s)  +  dFNis;  +  KNl)  drc(s-l)  +  K^M  df*(s-l)  +  dFM  (s) 


(2.46) 


in  which  the  pseudo  viscosity  load  vector, 


«/.  r»i  Jk  p«l 
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I?  (J  -j* 


d®(8"V 

(2.47) 


The  geometric  stiffness  matrix, 


°7a0  CNMQf&dA 


the  plastic  tangent  stiffness  matrix, 


Kir  -  hj£  +  12  jT  ^''W^dA 

( N ) 

and  the  nonlinear  pseudo  load  vector  dFM  (s)  contain  all  nonlinear 
terms  not  included  in  the  above.  It  should  be  noted  that  the  geometric 
and  plastic  stiffness  matrices  are  combined  with  the  responses  lagging 
behind  one  time  step  and  used  here  as  pseudo  loads. 

If  the  structure  is  fiber-reinforced,  then  the  suitable  transfor¬ 
mation  matrices  are  applied  to  (2.43).  The  derivation  of  governing 
equations  is  quite  complicated  for  asymmetric  loads.  Stresses,  strains, 
stiffness,  and  displacements  mdst  be  expanded  into  Fourier  series.  If 
loads  are  applied  statically, the  mass  term  in  (2.43)  drops  out.  For 

the  absence  of  viscosity  the  second  term  is  eliminated.  This  also 

(y) 

requires  the  pseudo  viscosity  load  dFN  (s)  to  be  dropped.  We  will  dis¬ 
cuss  these  special  topics  in  the  following  sections. 


2.5  FIBER-REINFORCED  SHELL 


In  general,  an  axisymmetric  shell  reinforced  with  fibers  would 
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normally  have  layers  of  angle  plys  with  each  layer  provided  with 
■f'crply  and  -cr  ply, a  being  the  wrap  angle.  Such  arrangements  tend 
to  keep  the  deformation  of  the  shell  still  axisymmetric  if  loads 
are  axiaymnetrically  applied.  A  typical  one  layer  element  is  shown 
in  Figure  A. 5.1. 

Orthotropic  properties  exist  in  the  directions  longitudinal  and 
tangential  to  a  fiber.  Coordinate  transformations  for  elasticity 
properties  are  then  necessary  to  conform  to  the  standard  generalized 
coordinate  system  of  a  shell.  A  similar  transformation  is  required 
for  the  plastic  stiffness  matrix.  These  transformations  are  elab¬ 
orated  in  Appendix  A. 5. 


2.6  ASYMMETRIC  DEFORMATION 


Let  g  represent  all  components  of  displacements,  strains,  stresses, 
stress  resultants,  etc.;  then,  we  have  the  Fourier  series  expansion 
of  g  for  asymmetric  deformation  of  the  form, 


g(«i8)  *  ^ 


cosjG  + 


2 

i.-i 


sinjQ 


(2.48) 


in  which  s  and  8  are  the  meridional  and  tangential  coordinates, respec¬ 
tively,  i  is  the  harmonic  number,  the  first  term  in  the  left  hand  side 
represents  harmonically  uncoupled  part,  and  the  second  and  third  terms 
denote  harmonically  coupled  parts  for  even  and  odd  functions,  respec¬ 
tively.  The  Fourier  series  representation  of  the  strain-displacement 


t 

x 


5 

\ 
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relatione  end  stiffness  matrices  is  given  in  Appendix  A. 6.  Equiva¬ 
lent  nodal  load  vectors  and  stresses  are  discussed  in  Appendices 
A. 7  and  A. 8,  respectively. 

2J  STATIC  LOADING 

2.7..  1  ELASTOPLASTIC  ANALYSIS 

There  are  in  general  two  methods  to  perform  an  elastoplastic 
analysis:  (l)the  tangent  stiffness  method  and  (2)  the  initial  stiffness 
method.  The  former  recalculates  or  updates  the  stiffness  matrix 
during  the  iterative  solution  cycles  whereas  the  latter  takes  ad¬ 
vantage  of  the  inverse  of  the  original  linear  elastic  stiffness 
being  kept  constant  on  which  all  subsequent  iterative  cycles  are 
based.  Obviously,  this  latter  approach  is  simpler  in  procedure 
but  requires  a  greater  number  of  cycles  for  convergence.  In  both 
cases,  the  load  is  applied  in  small  Increments.  The  present  study 
is  baaed  on  the  tangent  stiffness  matrix  approach.  Both  of  these 
methods  are  elaborated  in  Appendix  A. 9. 

2.7.2  VISCOELASTOPLASTIC  ANALYSIS 

Problems  of  creep  and  stress  relaxation  are  time  dependent 
phenomena.  The  procedure  for  analysis  is  identical  to  that  in  the 
elastoplastic  analysis  except  that  the  increments  in  loads  are  re¬ 
placed  by  the  Increments  in  time  and  the  total  load  is  applied  at 
each  time  increment  during  the  numerical  integration  of  the  govern¬ 
ing  equations. 

The  equations  of  equilibrium  are  solved  by  a  recurrence  for¬ 
mula  derived  from  a  suitable  difference  operator.  Such  a  formula  may 
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be  derived  by  assuming  a  linear  variation  of  either  displacements 
or  displacement  rates.  Explicit  forms  of  various  difference  opera¬ 
tors  are  given  in  Appendix  A. 10. 

For  each  time  increment  the  viscoelastic  analysis  is  performed 
initially;  then  the  plastic  matrices  are  calculated  from  the  visco¬ 
elastic  stresses.  A  standard  iterative  cycle  is  subsequently  re¬ 
peated  until  acceptable  convergence  is  achieved. 

2.8  DYNAMIC  LOADING 
2.8.1  ELASTOPLASTIC  ANALYSIS 

Effects  of  inQrtia  are  added  to  the  governing  equation  for  a 
dynamically  loaded  shell.  We  construct  the  mass  matrix  according 
to  the  formula  <2.39)  elaborated  in  Appendix  A. 11.  The  equations 
of  motion  thus  obtained  are  solved  once  again  by  a  suitable  diff¬ 
erence  operator..  In  general,  a  linear  or  constant  acceleration 
assumed  for  a  small  time  increment  is  used  for  deriving  a  recurrence 
formula  (Appendix  A. 10)'»  Various  schemes  of  numerical  integration 
have  been  investigated  and  are  available  in  the  literature  [2OJ. 

For  all  time -dependent  problems  as  indicated  in  the  previous 
section  the  Increments  in  loads  are  replaced  by  the  increments  in 
time  to  handle  nonlinear  plastic  behavior.  This  would  require  the 
total  dynamic  load  to  be  applied  at  each  time  increment  during  the 
numerical  integration  process  [2,3,4]. 


2.S.2  VISCOELASTOPLASTIC  ANALYSIS 


The  dytUufcic  analysis  of a.viscoelastoplastic  fiber-reinforced 
shell  is  our  ultimate  goal  in  the  present  study.  Because  of  this 
generality  we  will  summarize  and  itemize  the  analysis  procedure. 

(a)  During  the  first  time  increment,  the  generalized  displace¬ 
ments  are  calculated  from  the  recurrence  formula  with  initial 
and  boundary  conditions. 

(b)  Examine  the  current  viscoelastic  stresses  in  (a)  above  to 
determine  if  any  layer  through  thickness  of  any  element  has 
yielded.  If  so,  the  plastic  tangent  stiffness  coefficients 
for  that  portion  are  developed. 

(c)  The  plastic  tangent  stiffness  matrix,  if  non-zero  as  deter¬ 
mined  in  (b),  is  incorporated  into  the  recurrence  formula 
to  calculate  displacements  and  velocities  from  which  incre¬ 
mental  strains  and  strain  rates  can  be  found. 

(d)  The  incremental  equivalent  stress  is  calculated  and  compared 
with  that  of  the  previous  cycle. 

(e)  Steps  (b)  through  (d)  are  repeated  until  convergence  or  a 
certain  accuracy  is  obtained  for  the  incremental  equivalent 
yield  stress. 

(f)  For  a  yielded  element  (or  layer  of  the  element)  the  maximum 
equivalent  stress  which  was  originally  set  equal  to  the  input 
yield  stress  is  now  updated  by  adding  the  incremental  equiv¬ 
alent  yield  stress  to  account  for  strain  hardening.  For  the 
yielded  portion  the  anisotropic  parameters  of  plasticity 


must  be  updated. 

(g)  If  no  yielding  occurred  anywhere  in  the  shell  the  8tep3  (b) 
through  (f)  are  omitted. 

(h)  A  new  time  increment  is  initiated  and  the  above  steps  are  re¬ 
peated  until  desired  time  increments  have  been  completed. 
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SECTION  3 

THICK  SHELL  ANALYSIS 


3.1  GENERAL 

It  It  generally  known  [15]  chat  a  shell  with  Che  thickness-radius 
of  curvature  ratio  of  approximately  1/20  or  larger  is  regarded  as  a 
thick  shell  in  which  transverse  shears  become  significant.  Tradition¬ 
ally,  a  thick  shell  is  analyzed  similarly  to  a  thin  shell  with  trans¬ 
verse  shears  Incorporated.  Chung  and  Bandy  [6  ]  studied  finite 
element  applications  to  an  axisymmetric  thick  shell  considering  trans¬ 
verse  shears  in  the  planes  of  meridional  and  transverse  directions 
and  tangential  and  transverse  directions.  In  theiv  formulation^ rates 
of  these  transverse  shears  with  respect  to  meridional  coordinate  are 
added  to  the  bending  strains  of  respective  directions.  For  fiber- 
reinforced  shells  with  several  layers  of  angle  plys  through  the  thick¬ 
ness,  however,  the  three  dimensional  shell  theory  is  more  convenient. 

In  fact,  It  will  be  shown  in  Section  A  that  even  for  a  thin  shell  if 
layers  of  fiber  angle  plys  exist  through  the  thickness  direction  the  three 
dimensional  theory  appears  to  be  more  satisfactory,  particularly  when 
inelastic  yielding  is  considered. 

For  the  reesons  stated  above  the  present  study  will  employ  the 
three  dimensional  shell  theory  and  he  specialized  for  axisymmetric 
geometry  and  deformations. 
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3.2  INTERPOLATION  FUNCTION 

Although  the  three  dimensional  theory  is  used  only  radial  and 
axial  displacement  fields  may  be  considered  because  of  axisymmetric 
conditions,  thus  eliminating  shears  in  the  planes  of  axial  and  tan¬ 
gential  directions  and  radial  and  tangential  directions.  This  permits 
use  of  a  plane  element  subject  to  integration  around  the  circumference. 
SUcheplane  element  as  used  here  is  the  isoparametric  element  with  four 
corner  nodes  which  represents  a  linear  variation  of  axial  and  raditu. 
displacements  within  the  element.  This  element  was  originally  devel¬ 
oped  by  Ergatoudls  and  Zienkiewicz  [251.  Detailed  descriptions  of 
this  element  are  given  in  Appendix  A. 12. 

The  strain  components  considered  here  Include  normal  strains  in 
axial,  radial  and  tangential  directions  plus  shear  strains  in  the  plane 
of  axial  and  radial  directions  because  of  axisymmetry.  Development  of 
governing  equations  follows  closely  the  procedure  presented  in  Section 
2.  The  isoparametric  element  is  based  on  the  local  coordinates  with 
origin  at  the  center  of  the  element  and  four  corner  nodes  having  coor¬ 
dinate  values  of  1  and  -1,  This  permits  use  of  Gaussian  quadrature 
integration  in-,  the  plane  together  with  integration  around  the  circum¬ 
ference.  Consequently,  to  conform  to  global  coordinates  a  transforma¬ 
tion  by  means  of  Jacobian  matrix  is  needed  for  stresses  and  strains. 

Such  transformation  is  in  addition  to  fiber  transformation  as  elaborated 
in  Appendix  A. 5. 

The  generalized  displacement  field  ®t  (®|  ■  u,  ©#  *  v)  is  given 
by 


-  *1N®" 


(3.1) 
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where  i*l,2  and  Nal,2,...6  and  t1N  is  the  normalized  interpolation  func¬ 
tion  (isoparametric).  Alternatively,  8t  may  be  expressed  as 

Bi  -  ♦,«;  (3.2) 

where  r*l,2,3,4  and  tr  is  simply  the  interpolation  function  (isopara¬ 
metric)  given  by 

-  f(l-?)(l-'|) 

♦a  "  Jd+SHl-'n) 

+3  “  |(l+S)(l+T))  (i3> 

+4  -  ^(l-5)d+Tl) 

in  which  §  and  1)  are  Isoparametric  cootv  .«ates.  See  Appendix  A.  12  for 
derivation  of  (3.3).  The  difference  between  and  is  associated 
with  arrangements  of  thfe  nodal  displacement  vector. 

3.3  LINEAR  ELASTIC  CONSTITUTIVE  EQUATIONS 

Because  of  the  Isoparametric  coordinates  which  require  transfor¬ 
mation  to  the  global  cartesian  coordinates,  derivatives  of  any  variable 
with  respect  to  Isoparametric  coordinate  are  given  by 


where  the  Jacobian  J  is 
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Therefore, 


6r 

"1 

&z 

6? 

6r 

62 

W1 

b\ 

(3.5) 


(3.6) 


The  linear  elastic  stress-strain  relation  for  axisymmecric  solid 
is  given  by 

alJ  *  EU*A  Ykjfc  (3.7). 

where  t,4  •  1,2, 3,4. 

The  corresponding  strain-displacement  relations  for  small  dis¬ 
placement  using  (3.1)  through  (3.6)  may  be  expressed  as 
4 


YU  -  Yr  »  ~  -  2  n<A,n  +  ».  .5  +  C*  „Tl)ue 

Y»t  *  Y,  -  ~  -  2  fl(Arn  +  Brt?  +  CpnT|)vn 
4 

Y33  -  Y0  -J-  (1+B,..?+  C..1J  +  D,  5T1)u 


Y33  *  Y  - 


n 


(3.8a) 


(3.8b) 


(3.8c) 


Yis  -  Y„  - 


6u  6v  V  , 

+  Cn<Ara  +  B r„5  +  CrnT|)ut 


l>  T"  6z 


+  n(At8  +  B,,5  +  C,#f|)v#} 


(3.8d) 
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in  which  the  values  of  all  constants  are  defined  in  Appendix  A. 12. 


3.4  CONSTITUTIVE  EQUATIONS  FOR  FIBER-REINFORCED  VISCOELASTOPLASTIC 
SOLID (AXISWWETRIC  THICK  SHELL  APPROXIMATION) 

The  explicit  fora  of  the  thrae-dimensional. yield  function  given  in 
(2.6)  is  elaborated  in  Appendix  A. 2.  The  tensor  of  plastic  moduli 
fiiki  generally  referred  to  as  plastic  matrix  is  derived  in  Appen¬ 
dix  A. 2. 

Our  purpose  here  is  to  derive  first  of  all  the  orthotropic  visco- 
elastoplastic  stress'-strain  relation  of  the  type  (2.36)  specialized 
for  axisymmetric  three  dimensional  shell. A  detailed  procedure  achiev¬ 
ing  this  goal  has  already  been  presented  in  Section  2  for  a  thin  shelL 

The  final  fora  requires  change  of  indices  from  at ,p ,\,(i  to  l,j  ,X,A  so  that 

n 

dCT*J(s)  -  D1Jki  dykx(s)  +  C<A>d(qk.(s-l) 

Tm  1 


(r)  . 
+  B  dY. 


(r>  . 
+  C  dYk 


£(8)] 


JUJt 


dYkJe<8) 


(3.9) 


3.5  FINITE  ELEMENT  EQUATIONS  OF  MOTION 

In  the  case  of  a  thin  shel 1^ integration  through  thickness  was 
required.  Such  integration  here,  however,  is  combined  with  axial 
direction  over  the  area  of  iaoparametric  plane  element  by  means  of 
Gaussian  quadrature.  This  byt-isses  distinction  between  membrane  and 
bending  stresses.  Effects  of  f  .esses  are  combined  in  the  pre¬ 

sent  procedure.  The  final  form  of  equation  of  motion  following  the 
same  operations  described  for  a  thin  shell  in  Section  2,4  may  be 
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derived  in  the  form 


MM*d8"(a)  +  JN*d$(s) 


( • ) 

+  Kmm  d0(s) 


dFN(s) 


in  which 


<a)  (v)  (p) 

dFM(e)  -  dFN(s)  +  dFN(s)  +  dFN 


(r)  r  rCl  V  ‘1*4 

dFN(s)  -  [2nl  J  (  ?(r)  A  dc*»t j 

4l  tnl 


<r> 


+  Sty*  B  AHk£  AMl^r|j|d?dTl]d0H(8-l) 


<p>  r  rV1*, 

dFN  -  [2tt|  I  D 


J*4 


An i jAmu l  rjjjd^dll]d0(s-l) 


It  should  be  noted  that  A Hki  can  easily  be  Identified  in  the  expressions 
(3.8a)  through  (3.8d).  The  product  term  is  also  shown  in  expressions 
such  as  Fa#f  Gaa,  etc.  in  Appendix  A. 12. 


3.6  ANALYSIS  PROCEDURE 

Various  cases  of  analysis  including  elastoplastic  and  viscoelasto- 
plastic  analysts  under  static  and  dynamic  loadings  for  a  thin  shell  are 
described  in  Sections  2,7  and  2.8,  Other  than  three  dimensional  yield 
criteria  for  the  thick  shell  as  shown  in  Appendix  A.2fwe  follow  the 
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Identical  procedure  as  in  a  thin  shell.  Therefore,  no  further  com- 
ments  are  required  for  the  case  of  a  thick  shell. 
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SECTION  4 
APPLICATIONS 


4.1  GENERAL 

Bated  on  the  theory  and  procedure  described  earlier  various 
versions  of  computer  programs  have  been  written.  Detailed  descrip¬ 
tions  of  these  programs  are  given  In  Appendix  B  and  Its  suhappendices. 

Some  demonstrative  problems  are  solved  and  the  results  discussed 
In  the  following  subsections. 

4.2  GEOMETRIC  CONFIGURATIONS  AND  MATERIAL  PROPERTIES 

Three  types  of  geometry  are  considered  for  example  problems: 

(1)  cylindrical  shell,  (2)  combined  cylinder  and  spherical  dome,  and 
(3)  the  structure  same  as  (2)  with  a  portion  of  the  dome  section 
built  up  with  Increased  thickness. 

The  geometries  shown  in  Figures  1  and  2  are  analyzed  as  thin  shells 
whereas  those  shown  in  Figures  3,  4  and  5  are  treated  as  thick  shells. 

The  fiber  wrap  angle  a  Is  measured  from  the  circumferential  direc¬ 
tion  rather  than  from  the  longitudinal  axis.  Analytical  results  for 
isotropic  solids  and  fiber  angles  a  ■  0°,  50°,  90°  are  compared.  Deter¬ 
mination  of  lsotensoid  dome  surfaces  and  corresponding  fiber  angles 
based  on  netting  theory  is  Inadequate  if  bending  is  significant.  The 
analytical  means  of  determining  isotensoid  surfaces  by  bending  theory,  how 
ever,  appears  to  be  Impractical.  Furthermore,  wrapping  of  fibers  dic¬ 
tated  by  bending  analysis  in  general  may  not  be  possible.  For  these 


6.5 
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reason!  the  present  program  performs  stress  analysis  for  given  trial 
geometries  and  fiber  angles  to  determine  efficiency  of  performance  or 
load  carrying  capacity. 

In  our  examples  to  be  presented  in  the  following  subsections  two 
fiber  layers  are  used.  Each  layer  consists  of  angle  or  cross  plys 
wrapped  with  4Q t  and  -a  fibers  which  lead  to  axlsymmetrlc  deformations 
under  axlsymmetrlc  loads.  Material  properties  used  for  all  example 
problems  are  the  same  and  they  are  listed  in  Tables  1  and  2.  For  the 
analysis  as  isotropic  solid  the  material  properties  corresponding  to 
transverse  direction  are  used. 

4.3  STATIC  LOADING 
4.3.1  THIN  SHELL 

The  computer  programs  SP1  and  SVP1  are  used  for  analysis  of  the 
structures  as  a  thin  shell. 

The  elastoplastic  deformations  of  a  cylinder  (Geometry  Case  1, 
Figure  1)  subjected  to  internal  pressures  are  shown  in  Figure  6.  The 
only  experimental  data  available  for  comparison  is  plotted  also  in 
this  figure.  The  elastic  load  limit  for  analytical  solution  is  slightly 
lover  than  that  for  the  experimental  value  for  ar«50O.  It  is  interest¬ 
ing  to  note  that  for  ®"90®  which  represents  fibers  oriented  in  the 
axial  direction,  yielding  occurs  at  much  smaller  strain  and  no  strain- 
hardening  is  exhibited. 

The  viscoelastic  response  of  this  cylinder  with  c*»50°  subjected  to 
an  internal  pressure  of  250  pel  is  shown  in  Figure  7.  A  relaxation  time 
T(r)  *  .005  sec.  with  >•*'  1,2,3  and  the  integration  time  increment  At 


TABLE  1 

MATERIAL  CONSTANTS  FOR  THIN  SHELL 


Modulus  of  Elasticity 

El  ■  8.5  x  10P  psl 

E|  *  2.8  x  10P  psl 

Shear  Modulus 

Gli  •  1.3  x  10P  psl 

Poisson's  Ratio 

Vri  ■  .4 

M.i  ■  .131 

Plastic  Modulus 

e(pH  ■  8.5  x  103  psl 

E(p)T  ■  2.8  x  103  psl 
E(p)4no  *  3  x  10P  psl 

Plastic  Shear  Modulus 

G(p)tl  G(p)TT  -  1.3  x  1CP 

P8l 

Yield  Stress 

„ >l  ■  3.05  x  10P  pjl 
a(„),  -  3  x  103  psl 
^(0)46®  *  2  x  10*  psl 
°(o)TL  -  cr< o  > 1 1  ■  5  X  103 

psl 

Density 

.0388548  pci 

TABLE  2 

MATERIAL  CONSTANTS  FOR  THICK  SHELL 


Modulus  of  Elasticity 

Eu  -  8.5  x  lOPpsi 

E,  -  2.8  x  lOPpsi 

Shesr  Modulus 

Gjj  ■  1  x  1CP  psi 

Gli  "  1*3  x  MFpsi 

Poisson's  Ratio 

M.,  -  .262 

Vri  -  .4 

Plastic  Modulus 

e(j>)i  ■  8.5  x  10*  psl 

E(P)t  “  2.8  x  103  psi 

Plastic  shear  modulus 

G(p )l  “  1.3  x  103  psl 

G(p  >i  ■  1  x  103  psl 

Yield  Stress 

tf(o)i  "  3.05  x  10P  psi 

0(„ ),  *  3  x  103  psi 
<J(„)il  -  5  x  103  psi 

0(0  )tt  ■  2  x  103  psi 

Density 

.0388548  pci 
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■  .001  sec.  are  used  In  this  case.  It  Is  seen  that  the  viscoelastic 
displacement  gradually  increases  and  reaches  the  linear  elastic  dis¬ 
placement  at  approximately  .012  sec.  More  discussions  of  this  subject 
are  given  later  for  a  thick  shell. 

Effects  of  various  fiber  angles  for  a  dome-cylinder  are  shown  in 
Figure  8.  The  fiber  angles  for  the  dome  part  of  the  structure  are  arbi 
trarlly  designated  as  cy-S0,  15°,  25°,  35°  for  elements  1,2,  3  and  4, 
respectively.  With  this  dome  the  fiber  angles  for  the  cylinder  are 
taken  as  «■  0°,  50°,  and  90°  for  3  separate  analyses.  In  one  case 
Qf—  90°  is  used  throughout  the  structure.  Displaced  shapes  of  these 
4  different  cases  under  100  psi  are  compared  in  Figure  8.  For  a«0°  for 
the-  cylinder  the  displacement  is  largest  at  the  dome  but  smallest  at 
the  cylinder.  In  contrast  to  this,  <y»90°  everywhere  gives  the  least 
displacement  at  the  dome  but  considerably  larger  displacement  at  the 
cylinder.  For  cv«50°  for  the  cylinder  with  variable  fiber  angles  at 
the  dome  the  displacement  is  about  medium  throughout  the  structure. 
These  results  appear  to  be  quite  reasonable. 

Some  serious  difficulties  remain,  however,  in  the  elastoplastic 
analysis  of  a  fiber-reinforced  thin  shell.  These  difficulties  stem 
from  our  plane  stress  approximation  of  three  dimensional  shell,  which 
in  turn  affecr  a  unique  definition  of  anisotropic  yield  parameters  as 
proposed  in  Appendix  A. 2.  Specifically,  the  parameter  A1R  lacks  unique 
nees  as  the  fiber  rotates  between  0°  and  90°.  This  leads  Aj.;,  to  take 
on  large  values  which  consequently  cause  the  equivalent  stress  a  to 
assume  a  negative  value.  This  situation  arises  for  x  £  45°  in  static 
analysis  but  such  range  Increases  in  dynamic  analysis.  As  pointed  out 
in  earlier  sections,  the  plane  stress  approximation  for  a  thin  shell  is 


Figure  8. Displaced  Shapes  under  Internal  Pressure  of  100*  ppi', 
Variable  Fiber  Angles  at  Dome  and  Constant  Fiber  Angle 
al  Cylinder,  Static  Tinear  Clastic  Analysis  -  Geometry 
Case  2 
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simply  a  device  to  avoid  treating  a  ahell  as  a  three  dimensional  struc¬ 
ture.  it  Is  too  simple  an  approximation  to  afford  checking  the  yield 
phenomena  of  fiber- reinforced  shell  through  the  thickness. uniquely. 

Such  discrepancy  exists  also  In  dynamic  loading. 

4.3.2  THICK  SHELL 

The  computer  programs  SP2  and  SVP2  are  used  for  analysis  of  the 
structure  as  a  thick  shell. 

Figure  9  shows  plots  of  Internal  pressure  vs.  radial  displacement 
for  a  cylinder  (Geometry  Case  3t  Figure  3).  For  all  wrap  angles  other 
than  a*0°  the  elastic  limit  pressure  occurs  much  smaller  than  the  ex¬ 
perimental  value  for  50°.  The  shapes  of  the  curves  are  quite  sen¬ 
sitive  to  wrap  angles  as  noted  in  the  analysis  as  a  thin  shell.  The 
results  for  ar^O0  are  Identical  In  both  analyses  as  either  a  thin  shell 
or  as  a  thick  shell,  although  the  elastic  limit  for  <*-50°  In  the  thin 
shell  is  larger  than  that  for  ar*50°  In  the  thick  shell.  This  Is  at¬ 
tributed  to  most  probably  the  non-uniqueness  of  the  anisotropic  param¬ 
eter  Aia*  os  discussed  in  the  previous  Section. 

Displaced  shapes  of  cylinder  for  various  wrap  angles  subjected  to 
the  internal  pressure  of  2.27  psi  are  shown  in  Figure  10.  The  least 
displacement  in  both  radial  and  axial  directions  is  caused  for  a>0° 
whereas  or»90°  gives  the  largest  radial  displacement  but  very  small 
axial  displacement.  For  the  case  of  ar«50°  both  axial  and  radial  dis¬ 
placements  are  medium.  Of  course,  the  largest  displacements  occur  In 
isotropic  solid. 

Effects  of  integration  time  Increment  and  relaxation  time  are  sig¬ 
nificant  in  viscoelastic  analysis.  These  effects  are  even  more  critical 
in  viscoelastoplastic  analysis.  Figure  11  shows  effects  of  various 


Static  Internal  Pressure  vs.  Radial  Displacement  of  'Tode  1  (Elastoplastic) -  Geometry 
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a  -  90 
a  ■  0° 

a  -  50° 
Isotropic 


Displacement  Scale 
,  10“4  In. 


Figure  10.  Displaced  Shapes  of  Cylinder  (SP2)  for  Various  Wrap 
Angles  (90  -  a)  Subjected  to  Static  Internal  Pressure 
of  2.24  psl  -  Geometry  Case  J 


Fadial  isplacement  ar  Node 


Figure  11.  Effects  of  Kel.ix;ition  Time  on  Viscoelastic  Response 
of  Isotropic  ’ylinder  Subjected  to  Internal  Pressure 
of  2.275  psi  -  Oeometrv  Case  3- 
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values  of  relaxation  time  for  a  constant  Integration  time  Increment  of 
,001  sec.  tilth  an  isotropic  cylinder  subjected  to  Internal  pressure  of 
2.275  psi.  For  T(r)  *  .02  sec.  the  response  Is  extremely  erratic,  A 
alight  Improvement  is  obtained  for  T(r)  •  .001  sec.  A  stable  solution 
is  achieved  when  the  relaxation  time  is  reduced  to  .0001  sec.  at  which 
the  ratio  of  At  to  T(r)  is  10.  With  the  relaxation  time  being  the 
material  property  in  general  the  integration  step  site  must  be  adjusted 
for  a  given  problem.  Toward  this  end  Figure  12  shows  effects  of  var¬ 
ious  integration  time  increments  for  T(r)  “  .0001  sec.  For  the  fiber- 
reinforced  cylinder  of  or*  50°  and  isotropic  solid  subjected  to  Internal 
pressure  of  102.4  psi  integration  time  increments  of  .0002  sec.  and 
.001  sec.  are  used.  For  both  material  properties  At  ■  .0002  sec.  is 
unsatisfactory.  The  smaller  increment  of  At  ■  .001  sec.  produces  a 
stable  solution.  In  dynamic  analysis,  however,  it  is  well  known  that 
the  smaller  the  time  increment  the  more  stable  the  solution.  This  is 
an  interesting  contrast  associated  with  viscous  properties  of  the  mater¬ 
ial  which  tend  to  damp  out  a  motion.  Finally,  influence  of  ramp  load¬ 
ing  is  shown  in  Figure  13.  Once  again  a  proper  choice  of  At/T(r)  i® 

needed.  This  ratio  in  general  larger  than  10  appears  to  be  satisfac¬ 
tory,  although  yield  properties  may  require  .even  larger  ratio,  say  100 
or  larger. 

Stresses  are  of  fundamental  Importance  to  the  designer.  In  Figure 
15,  viscoelastic  stresses  in  directions  longitudinal  and  transverse  to 
the  fiber  are  plotted.  This  particular  analysis  corresponds  to  ar»50°, 
T(r)  ■  .0001  sec.,  and  At  ■  .001  sec.  in  Figure  12,  although  these 

stresses  are  calculated  for  the  element  1  which  contains  the  node  1  for 

which  Fig.  12  is  drawn.  It  is  interesting  to  note  that  the  viscoelastic 


—  . .  Stresses  In  Fiber  Direction,  oL 

- —  Stresses  In  Direction  Transverse 

to  Fiber,  o, 


a  *  50° 


1  2  3  A 

Time  (I0’asec.) 


Figure  15.  Viscoelastic  fiber  Stresses  in  Element  1  of 
Cylinder  Subjected  to  Internal  Pressure  of 
102.4  psl  -  Ceometry  Case  3 


stress  increases  until  t  *  1.5  x  10"3  sec.  and  decays  asymptotic 
to  the  linear  elastic  stress  at  approximately  3  x  10"3  sec.  Such  decay 
may  not  occur  for  hours,  days,  months  or  even  years  depending  on 
the  material  property  or  here  the  relaxation  time.  For  instance  con¬ 
crete  creeps  under  sustained  load  and  stresses  relax  over  the  period 
of  time.  The  relaxation  time  could  be  100  or  1000  times  larger  than 
used  in  this  analysis.  Such  large  value  then  would  require  corres¬ 
pondingly  large  integration  time  increments  to  maintain  the  ratio 
At/T(r)  for  stable  solution.  In  such  case  the  results  that  might  be 
plotted  in  Figure  15  would  simply  require  a  larger  time  scale. 

Results  of  analyses  of  a  cylinder-dome  (Geometry  Case  4)  and  one 
with  built-up  section  (Geometry  Case  5)  are  presented  in  Figure  14. 

For  the  cylinder-dome  it  is  shown  that  a*50°  contributes  to  consid¬ 
erably  less  displacement  than  the  isotropic  solid  does.  The  displace¬ 
ment  of  the  dome  with  built-up  section  is  very  small  at  the  dome 
portion  but  very  large  at  the  bottom  of  the  cylinder  portion  (Figure 
14b). 


4.4  DYNAMIC  LOADING 

The  computer  programs  DVP1  nnd  DVP2  are  used  for  analysis  of 
the  structure  as  a  thin  shell  and  a  thick  shell,  respectively. 

In  general,  displacement  of  a  structure  under  a  load  applied  dy¬ 
namically  is  60%  to  80%  more  than  that  under  the  same  magnitude  of 
loud  applied  statically.  Consequently,  dynamic  stresses  are  larger 
than  static  stresses,  and  design  of  the  structure  must  depend  on  ac¬ 
curate  or  realistic  limiting  stresses.  A  need  for  dynamic  analysis 
is,  therefore,  apparent  when  explosive  pressure  is  involved  such  as 


in  missile  structures. 


54 


The  linear  elastic  responses  of  both  isotropic  and  fiber-reinforced 
cylinders  (Geometry  Case  1)  subjected  to  Internal  dynamic  pressure  of 
250  psi  are  shown  Ln  Figure  16.  For  the  case  of  50°  considerably 
smaller  peak  response  is  observed.  It  is  also  seen  that  for  both  cases 
the  dynamic  peak  radial  displacements  are  approximately  80%  more  than 
the  static  radial  displacements. 

The  dynamic  viscoelastic  responses  for  the  cylinder  of  Geometry 
Case  3  subjected  to  dynamic  internal  pressure  of  34.1  psi  are  shown 
in  Figure  18.  Smaller  peak  responses  for  viscous  behavior  are  noted. 
Here  the  relaxation  time  of  10"*  sec.  with  Integration  time  increment 
of  lO'6  sec.  is  used.  This  gives  the  ratio  At/T(r)  of  1.  Such  ratio 
would  have  been  too  small  for  stable  solution  in  the  case  of  static 
loading.  For  dynamic  analysis,  however,  even  At/T(r)  "  1  is  very 
large  as  noticed  by  significant  damping,  particularly  for  Qf«50°.  It 
is  expected  that  the  ratio  At/T(r)  much  smaller  than  1  could  provide 
stable  solution. 

Figure  19  shows  the  dynamic  elastoplastic  and  vlscoelastoplastic 
responses  for  the  cylinder  subjected  to  internal  pressure  of  102,5  psi 
with  T(r)  ■  At  ■  10”6  sec.  Once  again,  effects  of  viscosity  contrl- 
buted  to  smaller  peak  responses.  For  the  case  of  a  •  50°  the  peak 
vlscoelastoplastic  response  is  considerably  smaller  thsn  that  for  the 
isotropic  solid,  but  such  difference  is  almost  absent  for  elastoplastic 
response. 

Dynamic  stresses  corresponding  to  these  peak  responses  are  al¬ 
most  twice  as  high  as  static  stresses.  Results  of  a  static  analysis, 
therefore,  would  lead  to  unsafe  design  when  such  load  is  to  be  applied 
dynamically  to  the  structure  in  service. 


Isotropic 

ur  -  r)0° 


At  •  10“*'  see. 


I’inure  lb.  T. inear  Kluatic  Response  of  Isotropic  and  Fiber 
Reinforced  Cylinder  Subjected  to  Internal 
I'reoaure  of  2r>0  pal  -  Geometry  Case  1 
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SECTION  5 
CONCLUSIONS 


Theoretical  formulations  and  computer  programs  for  the  analysis 
of  thin  and  thick  fiber- re in forced  arbitrary  axisymmetric  shells  sub¬ 
jected  to  static  and  dynamic  loadings  with  diversified  material  pro¬ 
perties  such  as  viscoelasticity,  elastoplasticity,  and  viscoelastoplas- 
tlclty.have  been  presented. 

Adequacy  of  anisotropic  egress -history-dependent  yield  parameters 
modified  from  Hill's  yield  criteria  and  internal  or  hidden  variables 
approach  for  viscous  behavior  has  been  demonstrated.  However,  the 
mlsotroplc  parameter  corresponding  to  a  plane  stress  in  the 

fiber-reinforced  thin  shell  is  not  unique  as  fibers  change  in  wrap 
angles.  It  is,  therefore,  necessary  to  analyze  the  structure  as  a 
thick  shell  in  which  a  three  dimensional  theory  is  utilized  if  yield¬ 
ing  of  fibers  is  to  be  involved. 

Dynamic  stresses  are  larger  than  static  stresses.  Although  de¬ 
sign  calculations  may  be  carried  out  with  static  equivalent  load  actual 
dynamic  responses  for  various  geometrical  and  material  properties  are 
too  complicated  to  be  guessed  at  from  so-called  dynamic  "load  factor". 

Vlacoelastic  stresses  are  large  initially  but  decay  as  time  elap¬ 
ses.  If  the  structure  is  under  sustained  load  creep  and  stress  relax¬ 
ation  may  be  important.  However,  selection  of  relaxation  time  must 
be  made  carefully  because  such  viscous  behavior  may  be  exhibited  for 
only  a  few  seconds  to  months  or  years. 
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APPENDIX  A.l 


STRAIN  DISPLACEMENT  RELATIONS  FOR  A  THIN  SHELL 


1 .  KINEMATICS 


Consider  the  location  of  a  point  In  the  undeformed  shell  defined 
by  the  position  vector 


r  *  *1  it 


where 


c,  =  (ft  »  f  »  ?  "  z) 


Thus, 


r  3  »,  +  zn  (1) 

where  is  the  position  vector  of  a  point  on  the  undeformed  middle  surface 
and  n  =*  n  (P1  ,  if)  is  a  unit  normal  vector  to  the  middle  surface.  The 
square  of  the  length  of  the  line  element  is  given  by 
(d£h)n  e  dr ’dr  =  ^pd^dF^-tfeadef  d!? 


where 

with 

Thus, 


dr  =  dr.  +  d  (zn)  *  dr  +  dzn  +  zdn 

*s.  m!y  *N. 


(2) 

(3) 


o',  p  *  1,  2  and  t,  j  3  1,  2,  3 


(ds,,  )a  ■  dr0  *dr<,  +  2dzn*drn  +  2zdn’dt^ 
+  za  dn*dn  +  dz;’ 


(4) 


Now, 


dr  *dc  =  r  dr^’t  QdP^  =  a  *a  d?f*dP  ■  a  d^dlj 

n  ut0  lnM  -  n  ,o  .•  'y  fl  ^  *  * 


3 


*n in  <y  "0'S' 


is  the  fundamental  form  of  the  middle  surface  where  a 


a  3 


the  first  funda¬ 


mental  tensor. 


O'  P  o'  B 

dn*dr„  *  n  *a„dP  dp  3  -b  dp  df  is  the  second  fundamental  form  of 

O  ru  Q  /-v/  ft 
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(ds)8  -  dR'dR  -  GtJ  dr*  d?J 

-  Gyp  d|a'd?P+  G^d^d«f+  (^d^d?3 

-  (A  +  zN  )  •  (A  +  zN  )df  d?8 
<*  -a  -  P  ~»P 

+  (A  +  zN  )•  Ndf  +dz8 
—  cv  ~,a  — 

■  %  •  2lB«>  +I*W  ^ 

+  (A  +  zN  )  •  N  dfdf  +  dz8 

SI  --,fY 


(7) 
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where 


%  -  U  ■  ■»„  +  *'% 

gm  '  ‘  <*«  +  *V  •  S 

033  *  1 

and*  Aap‘ V*p  •  ^  mb><*  •  ^  *5>>P 

B  «  N  •  A 

all  ~»a  -  P 

c  a  "  N„  *  N  a  "  <-B^  A.)  *  K  A) 

ap  ~,a  ~,P  a  ~X  p  -41 

Xu  1 

**  B  BP  A  -  B*  B_ 

a  p  AM-  a  PX 

The  difference  between  the  squared  lengths  of  the  line  ilements  on 
the  deformed  and  undeformed  surface  is: 


dsa  -  ds*  -  2Yj ,  d?|  d§J 


*  2Ya?  dfd?,P+  ZY^dfd?3  +  2Y33  dz* 

*  [(A  n  -a  )  -2z(B  -b  a)  +  zz  (C  fl-c  )]  dS^d^ 

aR  aR  ap  aP  a?  ap  s  * 

+  2(N  •  A  Q.  +  zN  •  N  ^  )  df  df 


so  that 


Y«s  ■  «<VV  -  22<%-  V  +  2S(c„e  • 

Y  =  N  •  A  +  zN  •  N 
Q?3  ~  -a  -  -,a 

Y33  *  0 


Denote  the  middle  surface  membrane  strain  as 


e  “  h  (A  .  -  a  a) 
aR  aP  ap 


A  ”  ®  (tb  +  u)  =  to  ,a  +  u  +u 

-a  -v  fit  ~  —  ,a  —  '  ~,a  -  -  a  ~,a 

A  =  A  «A  *a+a,y*u„  +  a*u  +  u  •  u 

ap  ~a  -  r  ap  -  -a  ~,P  ~R  -  ,a  .-fit  -,P 


(9) 


Thus, 


where 


e  .  -  k(*  *u  .  +  a  *u  +  u  *u  ) 

crp  --  a-  -~,3  —  P  — ,o t  —>a  —i  “ 


u*u^a.+tfJn 

~  ~  B 

u  *  (u^  a  )  +  if*  n  +  u3n 

-fit  -  p  fit  0-. 

=  a.  +  u^a.  +  u3  ,  n  +  t^n 


>ot—  3 


-  S  fit 


fit  -  -»Q ! 


in  which 


a  _  *  I"*1  a  +  Ja  n  *  r**  a  +  b  n 

~  PjCr  pry  -  |t  0a  —  Bar  ^ —  ja  By  ~ 


where  Christoffel  symbols, 

.41  I1 


„  -  a  „  *a  ,fa  =*  a  „  •  n  ■  a  *n  *■  b_ 

By  ~  g, ry  Pa  —  B,cr  —  ~  3  ~>u  3* 


Hence , 


.  ^  ,  lU  Px  ,  3.  i  3  3.^ 

u  =*  (u  +  1  u)a  +ubn  +  u  n-uba 

-,or  ,a  fhr  ~  I*  »B~  ,ot  *»  a  —  u 

=  (uM  +  rM  uS-u‘bM  )a  +  (uPb  +  u3  )  n 
,a  pa  a  -  n  ag  »a  - 

=  (u^  -  )a  +  (u3  +  u^b  )n 

ja  a  ~  u  »a  nfl  ~ 


where 


J  +  a* 

|a  ,a  3a 


and  the  "|"  represents  covariant  differentiation. 

Using  (13)  in  (10)  yields: 

°c 3  •  %f“-a  •  u"be  >?  m+  (u3»0  +  “M  W~ 

+  a  -[(of  -u3b*i  )a  +  (u3  +  u^b  )nl 

•  0  ja  a  ~  m  ,a  au  - 

+[(u^  -u3b^  )a  +  (u3  +  u^b  )n]  •[  (u^  -tf’b^a 

ja  a  -  u  ,a  au  ~  [3  3  -  u 

+  (u3,3  +  U"  b3u  >"  11 


(10) 


(ID 


(12) 


(13) 
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*  \  [u  ,  +  u  |  -Zt^b  „  +  (u  |  -t^b  )(u*f  tfb^  ) 

L  ale  e|a  ae  u|a  Ua  |p  B 

+  (tf*  +  A  )  (ua  +  uYb  Y  )i 

»a  aX  »B  3y 

Let  m  *  N  -  n  ,  N«m+n 
Denote  the  change  in  curvature  tensor  as 


-(B  -  b  )  ■  A  *  N  -a  *  n 


v  “  -d;*a  *  w 

*oB  ae  ae  -  a  -.g 


-a  -*e 


(a  +  u  )•  ni  -  u  •  bAa  , 

—  a,  ~,a  ~,B  ~»a  B  —  X 


Denoting  in  *  nr  a  +  mn  ,  then 

-  -  a 


since 


similarly, 


m  a  (m^,,  “  roh^  )a  +  (m  +  m^b  )n 

~>B  |P  B  -  U  »3  BU  - 


XaB  '  ma|B  ‘  ^  +  (“u|a  ’'‘V  )("|b  '  > 

+  (u!a  +  "\x  >(".B  +  ”Y  bBV  >  'C. ,, 


BY  '  B(u  i  ■«  >  ) 

K  n|a  vi 


ta- 

6°*  a  •  a  =  &  „  „  -  2n 

-a  ~P  *  aB  -  - 

i  aB 

n  s  $  £  a  *  a 
~  a  —  p 


n  =  %  A  -A 

*»  -  a  ~  p 


If  the  permutation  symbols  in  the  deformed  and  undeformed  surfaces 
are  the  same  (small  strain);  i.e.,  €°^  =  ,  then 


.  3  B.  ..XB  ^  /  Y  .  Y  3.  .3  ,  ,  n  ,  a 

f-u  -u  b  -  €  €  (ui  -b  u  )  (u  „  +  b  u  )}  a 

»a  aB  aY  X  X  »B  Bn  *" 
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+  f  u|a  -  ba  ^  +  ***  evA  <“|a  ‘b «  '  *£  “*>1  » 


We  define 


"a  *  ~u*a"u^bap  A*  ‘W  >  <“%  *  *>„/> 

"  ’  “fa  '  b“  "J  +  4  e<rf  6yx  (u|0  '  ba  "3>  ("|b  " > 

V  *(mH|a  •  ■*„  a  )(m|8  '  ■*?  >  +  <%  +  "W  V  ><%  +  > 


•  (Ma  ‘  mb\a  )bs  ’  <\|S  -  %>»;  1 


So  finally. 


Y  *e  +  zX.  +  z#A 

^  of  «P  aB 


The  shear  strain  is 


Y  -  N  •  A  +  zN  •  N 

op  -  -  a  -  -,a 

=  m  +m^u  .  +  (u3  +u^b  )(nri-l) 

a  Ufa  Ma  ,a  au 

+zmHn  .  -m^rnb  +zm(m  +ra^b  ) 

Ufa  iaa  ,a  au 

-zm^b  +  z(m  +  b  ) 

Ma  ,a  aM 


2<  SPECIAL  APPROXIMATIONS 


For  large  deflections  but  small  rotations 

,  B  3 
b  „  vr  «  u 


00 


.a 


/  a  ,  a  a »  3 

u  _  *  (uift  -  b  u  )a  +  u  n 

~»fl  |p  3  ~  a  ,3  ~ 
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Figure  A. 1.2:  Approximation  for  Small  Rotation 


a3 


i(u  ,  +  u  ,  -2b  u*3  +  u3  u*  ) 

a|g  p|a  ap  ,a  ,8 


m  ■  h  6°^  (2a  +  ud  n)  *  u3  n 

-a  »a  ~  »p  ~ 


from  which 


ah  «v  3  3  p  a  a 

*  k  (-26  0  an  u  =  -u  -  a  =  -u  a 

aP  —  »p  »p  >a 


m  =  0,  m  =  -u 

a  ,a 


Consequently: 


m 


rtJ 

1 

H 

a 

a 

-  xi' 

(r 

a  u 
n  a 

>a3 

— 

,a 

m3  - 

3 

a 

3 

.  a 

*  "U|a8 

a 

-  U  D 

,a  P 

n 

X  *  -(a  +  u*n)  •  (u*.  aY  +  u*  n) 

~a  »a~  |Y0'*  p  »Y 


-u 


Y  3  3 

*b„  u  v  u 


iae  wp  V  .a 


or,  neglecting  the  nonlinear  term^ 


X  = 


op  =  ’u|ap 


6* 


\ 


a  n  m  £  a 

fca@  a3 

ag  B  g03  1 
a 


€*  -  €*^  A 

aa 

6*op  m  g*cfi  i 


a  •  a 

-a  ~a 


s  B  n  ,  a  •  a; 
a3  ~  ~ 


n  •  a 
-  ~a 


a  ■ 


a 

a 


c  a 

* a a  -a 


,  n  *  a  ■  -€  a 
*  ~  -  U  ai>- 


a  ,  ^ 

ta  a  ,  a  *  n  -  -e  R  ay 
a3  ~  ~a  ~  ap  ~ 


Aa  ’ 


V  *  ?- 


-€  a 

Me  - 


a 


€^6  -  1 

aa 


Love-Kirchhoff  Assumptions 

(1)  points  which  lie  on  one  and  the  same  normal  to  the  undeformed 
middle  surface  also  lie  on  one  and  the  same  normal  to  the 
deformed  middle  surface. 

(2)  the  effect  of  the  normal  stress  on  surfaces  parallel  to  the 
middle  surface  may  he  neglected  in  the  stress-strain  relations. 

(3)  the  displacements  in  the  direction  of  the  normal  to  the  middle 
surface  are  approximately  equal  for  ail  points  on  the  same  normal. 

Assumptions  for  Love's  strain-energy  expressions 

(1)  the  shell  is  thin,  i.  e.  h/r  «  1 

* 

(2)  deflections  may  be  large  but  strains  are  small 

(3)  strain  energy  is  a  quadratic  function  of  the  strain  components 
*  -  |  eh^yijY kje 
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(4)  the  state  of  stress  Is  approximately  plane,  1.  e.,  the  effect 

of  transverse  shear  stresses  and  of  the  transverse  normal  stress, 

acting  on  surfaces  parallel  to  the  middle  surface,  may  be  ne> 

elected  In  the  strain  energy  density.  For  thin  shell,  terms 

with  z*  *  0  and  Ya3"  0.  For  shallow  shell,  if*  »  b  ftu^  »  0. 

»«  «P 

3.  CALCULATION  OF  FUNDAMENTAL  TENSORS 


a2 '  “  (alla22‘  ai2>  -Rjl+  (RJ  i 


ll  *  ^ 

|ax  Xa2 


(cob  9  i.  +  sin  Si,  -  R'i,) 
11  22  03 


‘“''ll..  2  <cos  9  i,  +  sin  9  i,  -  R'  1  ) 

r/TTr'2  -1  ~2  o  ~3 


o  o 


- "o  (cob  0  1+  sin  0  1,  -  R'  1  ) 

v/TTr'2  -1  o  ^ 


Sn  R'R" 
~  _  o  o 


s  ^7  *  _  2  (cos  0  1  +  sin  0  -  R'  1  )  +  — 2 

^3  /r  +  r'2  -1  0  ~3  /T 

RoRo'  R''(l  -  R'2) 

a  7rrr~2  (cos  9  £i  +  8in  9  io>  +  — - r~ 

/r+  r'z  -1  -2  /T7r'2 

O  o 

’  ^’/TTr'2  8in  9  -1  +  c°s  9  ^ 

o 

Rl2R"  R"(l  -  R'2)  R» 
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j 

i 

* 


R, 


b; 


iii 

’a 


(i  +  r:2)(/tt 


r:2) 


(l+R^2) 

~K 


3/2 


Since  a^  “  *  0 

parallels  are  lines 


R’2) 


-  R  /TT“R'2 
0  0 


we  say  that  the  meridians  and 
of  principal  curvature. 


*  Radius  of  cur¬ 
vature  of 
meridian  curve 
R  (x3) 


Figure  A. 1.4:  Principal  Curvatures 


R' 


o 


•  tan  <p 


AB  *  R  tan  cp  ■  R  R' 
o  T  oo 

BO.-'rJ  +  rZ  -  R„/l  +  RiJ  ■  -«2 


-  sign  signifies  R2  is  opposite  to  direction  n.  To  check  if 
a^p  and  R^  and  R2  define  is  valid  surface  we  use  Coadazi  equations 


V|y  "  Vie 


fa 

1  X  ,  22 


) 


1  /ST, 


Rj  v  "22"  '  Rj  ',1  ’  R^'ll'  '  Rj  '.2 
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_ ^ 


£i  " 
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2ix«a 


/alla22 
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-  a  •  n  b 
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/511  /a^ 


"2a  *  2.i  b?i 
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f*22  /a22 
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Similarly 


^*22  . 
**2m  *2  h 
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'11  ~1  ~,1 


/-  ^l  ®11 

/a.,  t,  *  - —  t,  -  - 


’ll  ~1  Rt  ~1 


R, 


4.  LINEAR  STRAIN-DISPtACEMENT  EQUATIONS  (NOVOZHILOV) 
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Let  Aj  a  .  A2-/T22  ,  R,  .  *2“-^ 


^  *  u^/a~ 


W* 


«3 


'/^cnaa03 


rfya  -  (r'u' +  r' u2)an 

Xs!  Qf|Ji  U  21  11 


M>  1 

u„a_  = 


8  ,3* 


•r^"/r;e  „  35s  1jt*  V 

M'M’ 


Baora  03 


W3$ 


a%  ^  V  w 


e  * 

11 


/a 


i  r  ,  .“l  “1  j/fl  — 1 

1 _ .III  ,  /—  <van  .  u  1 

~  /a  a  d§  11  -1  *  as.  ‘l/r-yr— 
n  il  L  li  1  Ai  1  1  /a.i  /flu' 


11  1111 


di/a~~  —2 

-t/sr;  t,  •  t  u 


ii  ^i  a?  i.i  . —  / - 

2  / a00  /a, ,a 


a  /alU-3 
11  (  Rt  )U 


22 v  "naii 

l  ;£7l  i  ^/®7i  —  2/577  -1  a,, 

1  du  1  ll  u,  +  y^-  - 11  u  11 


yi^"d5i  aU  d5l 


11  ft5‘  /*ii  /*u‘n  /an‘n 


VT,  -2  a 

+  /I7.  — 11  “  11 


_ _ _JL_  +  iL 

11  a52  /ro  /T-T-  /r~=~  R 


22  v  “11-11  /allail  "l 


u  *  u  ,  u  *  v,  u  *  w  '/*ii  *  **22  "  A 


-  m  1__  du  _ L_  w 

en  as1  ala2  3?2  Rj 


Similarly 


-  .  1_  5v  _J _2  ,  w_ 

C2  A2  a?2  AjA2  *1  R2 


„  tl  d_  .v_\  l  1_  /iL.\ 

u  Ai  as  a2  a2  a§  i 


-  i  av  ,  _i ii4  «_ 

22  a2  a^2  AjA2  a§j  r2 


\,b  "  ■  ‘  (“X,  ^'  ,a  -  b6<M«  '  “V* 

■  ‘  “|<*  •  “leV  •  u\ms  •  bs%k '  *V> 


3  .X  .X  .  n  .  ,  3 

-  'lot  -  "X|8b»  ■  “xb«j|S  -  V„|«  •  n5b1I~u 


i  a dw  u  .  _i _ ”*1,1  aw  v  . 

’  ax  a$x  \  as/  Rj  *  AtA2  a?,^A2  a^  '  rJ 


1  a ,J._  dw  v. .  _ 1_ ?  dw  u_ 

A2  ^  \  aS/  y  AtA2  aSj  Aj  aSj'  rx 


j. _ ,a2v  i  dw  _i_  ^2  aw  . 

a.a,  *a?,as0  ’  a*  a?,  ar,  *  a^  a?,  a?,/ 
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For  a  shell  of  revolution 
Let  Aj  =  1  A2  *  r 

bS 

*1  ■  -  sS 

«,  -  — — 

2  cos  Cp 


Novozhilov  strain  equation 

—  bu  bm 

e  =  —  -  w 

11  bs  bs 

_  1  ,bv 

e22  *  “  +  u  sin  <p  +  w  cos  rp) 


-  bv  v  1  bu 

e,0  *  t  -  ~*  sin  cp  +  ~  ■— 

12  bs  r  v  r  bG  ' 


"11 


v2  2 

b  w  b  m  bu  b'p  v 

-  (—7  +  »  — r  +  7T  T~  ) 
bs2  bs2  48  68 


<22 


1  f  1  b  W 
1  ,b2w 


COS  fD  bv . bw  .  •, 

“  be  +  sin  P(b7  +  u  £  )} 


12 


bw  Dq>  bu 

be  "  bs  bo 


"  7  ^tsbO  ’  f  sin  <p  -  rr  ~  +  cos  ff>  rr  - 


bv  cos  <P 


bs 


sin  w  } 


(16) 
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APPENDIX  A. 2 
YIELD  CRITERIA 


I.  ISTROPIC  YIELD  FUNCTION 

The  von  Mises  yield  criterion,  one  of  the  most  widely  used,  is 
discussed  herein.  Von  Mises  suggested  that  yielding  occurred  when  the 
second  deviatoric  stress  invariant  J  reached  a  constant  value  o  called 
equivalent  stress  such  that 

f  “  a3  =  3J  =»  |  T*ot  (1) 

where  f  is  the  plastic  potential  function,  and  T00t  is  the  octahedral 
stress  defined  as 


oc  t 


»aa  V  +  (0;,3-o u  )9  + 


*  a  a 

6(aj8+o>;t+<J3  ^  ) 

(2) 


Together  with  Drucker's  postulate  that  the  yield  surface  is  con¬ 
vex  and  the  plastic  strain  rate  vector  is  normal  to  the  yield  surface 
the  Prand tie- Reus s  flow  rule  for  Isotropic  hardening  material  as  de¬ 
fined  in  (2.8),  Section  2twill  he  employed  in  the  present  study. 

Following  the  procedure  outlined  in  the  expressions  (2.8)  through 
(2.19),  we  obtain  the  explicit  physical  components  of  the  tensor  of 
plastic  moduli. 


*  -D  Z  E1  D 

~  =  E<V>  +  ZT  DZ 


O) 


where  E  is  the  standard  elastic  matrix  and 
,  •  , 

-^91  3°33  3°lz  3^83  3°3X 

~T  =  *•  2<T  2<T  IF  7  ~o 

in  which  ,  a99,  and  <j3:j  are  deviatoric  stresses. 


(4) 


2.  ANISOTROPIC  YIELD  FUNCTION 

2.1  THREE  DIMENSIONAL  SOLID 

Developments  of  anisotropic  yield  criteria  have  r.voived  for  the 
last  two  decades.  The  directional  anisotropy  was  studied  by  Hill  [8] 
using  the  von  Mises  yield  condition  and  its  flow  rules.  He  chose 
yield  functions  such  that  the  von  Mises  criterion  for  isotropic  solid 
is  restored  when  anisotropy  vanishes.  Such  yield  functions  were  used 
also  by  Tsai  [22],  Fulton  [7  ]  and  others.  The  formulation,  however, 
is  not  convenient  for  numerical  scep-oy-step  computation.  Hu  [loJ 
extended  the  Hill  theory  and  Introduced  equivalent  stress-strain  re¬ 
lationships  in  appropriate  form,  using  simple  uniaxial  and  shear 
stress-strain  tests  on  coupons  cut  in  the  directions  of  the  ortho- 
tropy.  The  anisotropic  parameters  were  held  constant.  However, 
experimental  evidences  show  that  anisotropic  parameters  fnr  strain- 
hardening  material  depend  on  the  state  of  stress  and  should  be  up¬ 
dated  as  the  stress  level  changes.  Such  procedure  was  used  by  Jensen 
[ll]  and  Whang  [23].  In  the  present  study  the  yield  criterion  of  Hill 
as  extended  by  Hu  and  further  by  Jensen  and  Whang  will  be  specialized 
for  fiber-reinforced  shell  structures.  Additional  discussions  con¬ 
cerning  coordinate  transxotmatious  ^,:e  given  in  Appendix  A. 5. 
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For  the  case  of  anisotropy  the  yield  function  (2.6)  is  given 
explicitly, 

f  =  0s  =  —  f Aj *  (on  “Oga  )*  +  A«3  (oag-033  )*  +  Aaj  (ct33 -o^  j  )® 

+  6(A440X3  +  Aa>°*a  +  Ajr,af1)}  (5) 

where  AlB,  etc.  arc  anisotropic  parameters. 

Setting 

An  *  Ala  +  Aax 
Ass  "  Axa  +  Aaa 
Asa  *  Aaa  +  Aai 

f  *  a  =  J  {  AuO^  +  Aaaoaa  +  Aas^aa  *  2(Aiaanaaa  +  Aa30a#cr33 

g  J  2 

+  Aii^aa^xx)  ^CA+i0!*  ■*"  AiB^sa  Aon^.n)^  (6) 

Differentiating  (6), 


d,JXX 

"  A.'l  l(,33  ) 

=  2o 

(au°j  1 

”  Ato°aa 

doaa 

+  To” 

(■Ais°ii 

+  Aaaoaa 

1  "  Aa3CT33  ) 

do33 

+  -25- 

("A3XaXX 

‘  a»30«» 

A33°33  ) 

dolt 

da#3 

do3i 

+  (3^4^11)  ^  (3A B5^®3)  “  (3^6^31  )  (7) 

CTO  O 
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It  is  seen  that  (6)  reduces  to  the  isotropic  yield  function  (1)  if 
An  *  Aa0  *  A33  *  2  and  all  other  anisotropic  parameters  are  equal 


to  1. 

The  anisotropic  parameters  initially  depend  on  the  initial  yield 
stresses  in  various  directions.  Let  CT(0)  be  the  equivalent  initial 
yield  stress  in  the  1-1  direction  and  all  other  stress  components  be 
equal  to  zero,  and  set 


'H  0)  *  fT(  o)n 


which  gives  from  (5) 


«*<o)  *  '>(A(oHj>  +  A(  o)3i)a(o)u 


A(  0  )ia  +  A(o)3X  *  2 


Similarly  for  the  2-2  direction, 


°SP  =  °(o)SS 


-  2  1/  .  a 

°(o)  *  ^vAjoixa  +  A(  o  )g;i  )<)(  0  )gg 


A;  0)1  a  +  A(o)*..,  *  2 


CT(  0 )aa 


<?33  *  <7(0)33 


-*1  2 

a(0)  “  7  (A(o)23  +  A( o )3l )a( 0)33 


A(o)a3  +A(0)ax  *  2  57^7 


Solving  for  A(0>i2,  A(0)33»  and  A(0)3l  from  (9),  (10),  and  (11), 
get 

ja 

A(o)xs  s  1  +  °(o) 

_a 

A(  o)aa  =  -1  +  0(  o) 

J  I 

A(0)at  =  1  +  a0 


°(  0  >33  "  °( 3 )as 


°fo>33Crfo)aa 


a  s 

[<3(  0 >33  +  o)*al 

a  a 

a(  0>33°( 0)08 


°< o)aa  "  CT( o)a3 
5  i 

<3<  0>33a(  0>38 


\ 


Similarly,  we  determine  A(0>44»  A(0)5!5i  and  A(o)e(;  as 


Let 


A(  0>44  *  "3 


0(0)  \a 


a(  o)iaj 


/  - 


1  °( o) 

A(0,s"  '  3  ^TuT 


A(  o)> 


,J(o) 


M  1  ’  (  o )  a: 


A(0)ix  *  A(0)la  +  A(0)  31  *  2 


A(o>aa  *  A(0)ia  +  A(0)ta  “  21 


P(o) 


A(o>33  *  A(0)33  +  A(o)31  *  2 


a( o)aa 
?(o) 


f-M 
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These  initial  parameters  do  not  remain  constant  as  the  material 
undergoes  strain-hardening.  The  subsequent  anisotropic  parameters 
should  depend  on  the  current  state  of  stress,  equivalent  stress,  and 
bilinear  plastic  moduli  for  various  directions.  Wc  make  use  of  the 
fact  that  plastic  work  performed  in  t lie  current  stress  space  and 
equivalent  stress  space  for  a  given  direction  must  be  equal.  For 
the  l-l  direction, 

<p)  (*  (p)  r  _(p) 

W  *=  I  Ulldvll  *  I  JdV  , 

(p )  f  (p)  f  -  _(p> 

W  =  J  atidY„  -J  cd7  , 

etc. 


Figure  A.2.1:  Plastic  Work 


From  the  figure  above  it  Is  easily  seen  that 


eT’»VxT’(0<  •>l*‘  '  a(0)u)  +  2E(p)U  ‘  a<oHl) 


<J(o)a*  i  s 

E(  ,),t  (<7(,)”  '  a(0^«>  +  2E(p)„  ‘ 


9  E(p)»*  p  3  n  .  a 

a(«)8B  *  E(  —  l°(  On"  <’(  o)n  J  +  a(  o)»a 


Therefore,  the  subsequent  anisotropic  parameters  are,  using 


I  s  \* 

.  2  - 

°(  ■  )  1 1 


Aaa  *  2 


a(a)ia/  E<P>aa  (5a  .  o?o)u)  +  of0)ta 


E(  p)tl 


Similarly, 


E< p)aa  _o  a 

(°  ‘  o )  1 1 )  +  °(  o)aa 


G<  p)ia  _#  a  a 

3  I  _  °<  o)u  >  +  °<  oHal 

p)ll 


p )  l.i  r>  2 

3  i"F7~Z - (°  -  °(o)u>  +  °(o)ia  1 
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o 

A,  s>  .  - .  '  ■  ■  ■—  '  —  —  ■ 

G(  p)aa  *  a  a 

3C  -  («  *  °<0)il)  +  «(0)»3]  d**) 

■M  p)u 

where  E( p)  end  G( p)  arc  the  plastic  modulus  and  plastic  shear  modulus 
in  bilinear  stress-strain  curves. 

2.2  PLANE  STRESS 

For  the  case  ot  a  shell  subjected  to  the  state  of  plane  stress, 
the  yield  function  is  written  in  the  form, 

9  Au  r>  A**  s  2 

f  *  °  m  ~2~  all  ‘  Aii°Uaa»  +  —  °s»  +  3A**fIxa  (15) 

Here  we  can  only  provide  two  tensile  tests  and  one  shear  test 
which  arc  not  sufficient  to  determine  4  parameters  of  anisotropy.  To 
settle  this  problem  the  yield  stress  in  the  direction  of  thickness  of 
plate  or  shell  may  be  assumed  to  be  the  same  as  that  of  1-1  or  2-2 
direction  [llj.  This  assumption  was  rejected  by  Whang  [23]  who  sug¬ 
gested  additional  tensile  test  at  some  angle  from  the  axis  of  ortho- 
tropy.  In  the  present  study,  the  later  procedure  is  followed.  If  this 
tensile  test  is  performed  ot  an  angle  0  from  one  of  the  axes  of  ortho* 
tropy,  then  (15)  becomes 

_a  Au  a  a 

°  *  “  •»„  c°*\  •  Aia'>„  cos2()  sin  2Q 

Afa  2 

+  — 2 —  8in4((  +  3a44ou  cos  2q  sin20  (16) 

o 

For  convenience  the  tensile  test  may  be  performed  at  0  ■  45  which 


will  simplify  (16)  such  that  for  initial  yielding, 

j3  A(0)n  A(0)i8  A(0)#f  3 

°(o)  *  (7<0)4S°  l  — g - 5 +  8 -  +  4 

or 

A(o)u  A(0)ii 

A(0)i2  *  — 2 -  +  — 2 - +  "  * 

'  / 

where  A(o)u>  A(0)as?,  A(0)44  are  the  same  as  in  the  three  dimensional 
case.  The  subsequent  parameter  An  for  strain-hardening  is 


0(0 


a( 0 )4fi 


(i7: 


AU  Ata 

Aia  *  ~Y~  +  "2*  +  3A*4 


E(oHb° 
E(  p)n 


4o(0) 

a  a  a 

(o  ‘  J(o)u>  "  °(o)46° 


(18) 


where  An,  Asa>  and  A44  are  the  same  as  in  the  3  dimensional  case. 


2.3  ANISOTROPIC  PLASTICITY  MATRIX 

With  all  anisotropic  parameters  defined  it  is  a  simple  matter  to 
apply  the  associated  flow  rule  to  obtain  the  anisotropic  plasticity 
matrix  in  the  form 


*  -D  z  ZT  D 
D  ~ 


D(,)u  +  zT  D  z 


(19 


where  the  components  of  Z  are  given  by 
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3.  ANALYSIS 

The  linear  elastic  analysis  is  first  performed  to  calculate  dis¬ 
placements  and  stresses.  If  the  stress  in  any  element  exceeded  the 
limiting  yield  stresses  then  the  plastic  matrix  as  derived  above  will 
be  calculated. 

For  the  case  of  a  thin  shell  Integration  of  plastic  stiffness 
matrix  must  he  performed  with  each  layer  integrated  one  at  a  time  and 
summed  through  the  thickness.  Necessary  equations  for  this  process 
are  derived  in  Appendix  A. 9. 

Details  of  elastoplastic  analysis  are  also  given  in  Appendix  A. 9. 


APPENDIX  A. 3 


DERIVATION  OF  INTERNAL  (IHDDDI)  VARIABLES 
Consider  Che  internal  variable  * <j| } 

"I  eXPf^o)  ^*(T)dT 

where  V14(T)  nay  be  considered  to  vary  linearly  within  the  email 
time  interval  At, 

v,  J<T>  *  V./t-At)  +  T-^tAt)r»ti(t)  -  Vu(t-At)] 


Substituting  (2)  in  (i). 


<t>(t)mf  eXPt^T~j  v,i(T)dT  +  f  exp(-~~)  VtJ(T)di 
J0  *  •'t+At  r 

/-At  \  (r)  fL  /-(c-tX  . 

•  ™pfe7,  «-.“•»•> +l  cxp  (t~t)  v,i<,>dT 

x  7  •'  t-At  %  7 

'  ‘"p(€))",:‘(“At>  +J  ",p 


At-t+r  . 

+  [vu(t)  -  VIJ(t-At)J)dt 
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+  T 


(r ) 


•  “"(if:,)  >  'V'-a') 

<t)J  -  A  ,u(t-it) 


+ 


( r ) ' 

B  V, + 


(  r  ) . 
C  V 


1 J 


(t) 


or 


At)) 
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APPENDIX  A.  4 

INTERPOLATION  FUNCTIONS  FOR  A XI SYMMETRIC 
THIN  SHELL  ELEMENT 

Let  us  consider  an  arbitrary  axis ymme trie  shell  as  depicted 
in  Figure  A.4.1.  The  nodal  circles  passing  through  the  nodes  normal 
to  the  surface  are  boundaries  of  each  element.  The  deformation 
state  of  each  discrete  element  (Figure  A.4.2)  is  described  by  four 
midsurface  displacements  (Ref.  [24]). 

:i)  meridional  translation  ^ 

(2)  circumferential  translation  0* 

(3)  transverse  translation  ^ 

(4)  meridional  rotation  04 

The  generalized  displacements  at  nodes  p  and  p  +  1  for  the  ele¬ 
ment  with  meridional  length  l  are: 


at  s  « 

0 

at  s  •  l 

0i.p 

■  u(p) 

®l,pfl 

1 

C 

®s.p 

*  v(o) 

®e  p*l 
> 

> 

B 

9:, ,  p 

■  w(0) 

®3,P*l 

-  w(£) 

04  »  p 

-&♦ 

«jS*> 

bs'c 

®4  ,P'  l 

bw  &cp 

where  fp  is  the  angle  between  the  z  axis  and  the  line  tangent  to  the 
meridional  surface,  s  is  an  arbitrary  distance  from  the  node  p  along 


the  meridian. 
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The  mldsurfece  meridional  slope  cp  is  assumed  to  vary  as 
<p(s)  *  a0  +a^  a  +  8*8* 

Assuming  that  the  curved  shell  element  may  be  approximated  as  a  seg¬ 
ment  of  a  circle,  the  coefficients  a0,  ax,  and  a„  are  determined  as 


a0  -  cp, 

-  j(&Pi  -  4«,  -  2cp,fX) 

as  -  7e(*Pi  +  3ePp  ’  6«pL) 


The  meridional  variation  of  midsurface  displacements,  u,  v,  and 
w  may  be  assumed  to  be  of  the  form 

wx(«)  ■  u(s)  -  C;  +  CgS  (la) 

®»<«)  "  v(s)  -  C3  +  C4s  (lb) 

%(»)  ■  w(s)  ■  Cr,  +  Q;  s  +  C7sb  +  CflSB  (lc) 

The  meridional  rotation  is,  then, 

04 (s)  -  +  u  7?^  -  Co+2C7s+3Cn8;?+(Cx+C2s)rf,/(ld) 

Q8  0® 

Writing  these  four  equations  at  p(s«0)  and  p-H(s*j2)  we  obtain  8 
equations  from  which  we  can  solve  for  the  8  constants,  through 
Ce.  Substituting  these  constants  into  (la)  through  (Id)  we  have 


9i  *  s^q,,  eN 

(2) 

or 

®»  *  *1N 

(3) 

In  matrix  form , 

0«sQw«tw 

*  **  •  <•»  N. 

(4) 
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* 

where 


£  -  S  Q 

(4  x  8)  ('  x  8)  (8  x  8) 


(5) 


hi 

*  *M  *  1  -  f  . 

*18  "  *»0  “ 

s_ 

7 

'hi 

B  .  2CPqS8 

-  -SCPo  +  - 

-  t 
i  * 

»  i  -  ■ 

3s*  2s3  . 

g*  +  jJ3  »  *»4 

=  s  ■ 

2sa  A  s3 

•t  +  f 

*38 

i  i9  * 

.  Isi 

*>■'  *  t 

i 8  » 

-s*  s3 

> 

i 

&|  ** 
srl 

+ 

>° 

a 

68 

"  “  j e8 

+  .  *.4  •  1 

4s 
"  l 

+  Isl 

+  J?8 

*46 

>  *»  *  f1  - 

6s? 

F"  * 

u 

00 

* 

! 

ell 

other  *lN  -  0  , 

■*>'  *i“l- 

etc. 

Thus,  the  strain-displacement  relation  may  be  written  in  the  form, 


Y<*p  "  eap 

The  engineering  strain  vector  in  matrix  form  considering  only  linear  terms 
can  be  obtained  by  substituting  the  Interpolation  functions  into  the  strain- 
displacement  relations, 


Yi  -  ^(Cx  +  C8s)  -  (Cfe  +  Qjs  +  O7 8*  + 

b8 

+  C{  «*  +  Gb 8  +  c7s"  +  CHS8) 

+  (cx  +  0,8)  |3fS-  +  5_  (Ct  +  Cas)  } 

bs 

Writing  similarly  for  Ya  and  Via*  performing  differentiation,  and  sub- 
stituting  values  of  constants,  we  obtain 


X"£+C&",£H..9®+zwq.0 

e  -  GW  58 

Ci-z»S® 

where 


1  0  0  0  0  0 
G  -  0  1  0  0  0  0 

0  0  1  0  0  0 


(6) 


The  explicit  form  of  W  will  be  shown  in  Appendix  A. 6. 


APPENDIX  A. 5 


* 


TRANSFORMED  ELASTICITY  AND  PLASTICITY 
MATRICES  OF  FIBER-REINFORCED  AXISTMMETRIC  SHELL 


1.  GENERAL 

Two  types  of  element,  plane  stress  element  of  axisymmetrlc  thin  shell 
and„l8oparamatric  element:  of  axisymmetrlc  thick  shell  will  be  discussed. 
Shear  stresses  in  the  plane  of  meridional  and  tangential  directions 
are  included  in  the  axisymmetrlc  plane  stress  shell  element  whereas 
shear  stresses  in  the  plane  of  axial  and  radial  directions  are  incor¬ 
porated  in  the  isopcrametric  sdlid  element.  Coordinate  transformations 
between  the  local  and  global  systems  differ  depending  on  the  choice  of 
either  plane  stress  element  to  be  used  for  a  thin  sheil  or  isoparametric 
element  for  a  thick  shell.  Since  the  global  coordinates  for  the  thick 
shell  consist  of  axial  and  radial  directions  additional  transforma¬ 
tion  for  an  inclined  element  is  required. 


2.  PLANE  STRESS  ELEMENT 

It  can  be  thown  that  the  transformation  between  the  local  and 
global  components  of  strain  for  a  thin  shell  is  related  by 


cos*  a 

-2sinofcoacs' 

-sin*0 


cos<y  sin  a 
cosaof-8ln^ 


-coact  si  no 


os  pot 


(1) 


^WHWTiWnwiin. . .  ■....  r  .  - 


’’.'^wgja^jy-yiwwq  jy.wpgra^.ww  jwy 


■wTfv -’'.^7'' 


swsw 
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The  above  notations  are  shown  in  Figure  A. 5.1.  In  matrix  form  (1) 
may  be  written, 


(t)  (a) 

Y  =  T  Y  n) 

where  (i)  -2nd  (a)  refer  to  "local"  and  "global",  respectively.  It  is 
s  simple  matter  to  show  that  the  global  elasticity  matrix  is  in  the 
form. 


(a )  .  U) 

D  -  Tt  D  T  (3) 


U>  <i)  (  l  ) 

in  which  D  •  D(#  )  +  D(  p)  with  (•)  and  (p)  denoting  "elastic"  and 

(a) 

"plastic".  The  components  of  D  are 


(a) 

-  c"(caD,x  +  s3Dla)  +  4sacaG  +  8a(ca»xt  +  s8^,) 

(a)  „  . 

D;a  *  sc(caDu  •:  eaD13)  +  (c8-88)(-2scG)  -  cs(c*Dlt  +  eaD#a) 

(a)  .  . 

DX3  *  s,(caDu  +  8*Dla)  -  4s8cBG  +  c*(cBDlB  +  S*bts) 

(a) 

Dax  *  c38(Da-DV!,)  +  <ca-s8)G(-2sc)  +  CS3(D18-I%t) 

(a) 

D„  *  casa(Dll-Dl,)  +  (ca-8* )G  -  c8s8(Dlt-D„) 

(®) 

ftsa  *  c83(DX|-Dl3)  +  2sc(c8-s8)G  +  C3 s(Dia  ’D*a) 

(a) 

Ebx  *  c8(s8Dxj  +  c8Bxa)  -  4s-c*G  +  #B(s8Dxa  +  c*Daa) 

(a  > 

l>3t  <=  sc(s8Du  +  caDxa)  +  2scG(c8-s#)  -  sc(svD|9  +  c8Dat) 

(a) 

Eba  “  89(8aDu  +  c8^,)  +  4s*caG  +  c,(s,Dll  +  c8!^*) 


where  c 


cos  or  and  s  *  sine*  .  For  axlsymmetrlc  deformations  we 


%  r-  Meridional  displacement 
0n  :  Tangential  displacement 
f«i,  -  Transverse  displacement 
^  Meridional  rotation 


Figure  A. 5.1:  Fiber-Reinforced  Axlsyanetric  Thin  Shell 
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a  d 

require  that  stresses  carried  by  a  fibers  a  p  be  the  same  as  stresses 

(a) 

rv8 

carried  -o'  fibers  .  This  is  satisfied  by  imposing  the  average 

stress  to  assume  the  form 


+o*P  ) 
2'  <»>  <-»)' 


However,  the  strain  in  each  ply  is  equal  to  the  Imposed  strain.  This 
(o)  (a)  (a)  (q) 


causes  Dxo 

(a) 
of  D  . 


D»j 


DS:, 


0  and  we  obtain  the  following  form 


(a) 


•  (a) 

0 

(a) 

Pl3 


0 

(a  ) 
Oat 


(a) 

I>13 

0 

(a) 

^33  • 


(4) 


which  is  the  globally  transformed  fiber  reinforced  elastoplasticity 
matrix.  If  the  plasticity  matrix  is  used  as  pseudo  load  the  above 
operations  are  performed  separately  without  adding  together. 

3.  PLANE  STRAIN  ELEMENT 

Referring  to  Fig,A.5.2  the  coordinate  transformation  of  the  elas¬ 
ticity  and  plasticity  matrix  for  a  plane  strain  axisymmetrlc  solid  ele¬ 
ment  for  a  thick  shell  without  bending  may  be  accomplished  in  the  same 
manner  as  in  the  plane  stress  case  if  such  transformation  is  limited  to 
a  plane(vertical  cylinder).  The  global  elastoplasticity  matrix  cart  also 
be  given  equivalently  by 


(6)  .  <l) 

D  -  Tt  D  T 

4x4  4x6  6x6  6x4 


(5) 
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(t) 

in  Which  £  is  Che  sum  of  regular  3-dimensional  elasticity  and  plas¬ 
ticity  matrices  and  T  is  related  by 


K' 

l 

0 

0 

0 

vT 

0 

c* 

8* 

0 

Yl 

0 

8* 

C* 

0 

Ypt 

m 

0 

0 

0 

c 

Ypl 

0 

0 

0 

-8 

Ytu 

0 

-2cs 

2C8 

0 

— 

where  c  *  c.osa 

or 


s  ■  sir** 

<e>  (  .s 

y  »iyU) 


(6) 


If  the  element  is  inclined  an  angle  ofrp  from  the  z  axis 
(FigA.5.3),  then  T  must  be  modified  as  follows: 


V 

10  0 

c  0  -s 

‘•I1" 

M 

0  8  c 

0  10 

>1* 

? 

0  c  -s 

s  0  c 

1  a 

where 

F  c  0  -5 


«8 

S 

cc 

-  88 

c 

-sc 

and 
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c  ■  cot  {p  ,  a  ■  sin  (D 

Tj»  -  r  ,  71*  -  0  ,  if  -  z 


«a  defined  in  Figure  A. 5.2. 

The  modified  T  then  assumes  the  form, 


? 

*• 

0 

-S3 

c*i» 

c*c* 

8* 

c*8c 

a*sB 

a*5* 

c* 

s»3c 

2c3c 

-2BcC 

0 

B*c-3*c 

-2cei 

2alc 

0 

-^S+38S 

-2ci#a 

«2c^cs 

28C 

>2ccsB 

Here  T  replaces  T  in  (5)  for  Inclined  element. 


Ftl'iire  A.r).|:  Coordinate  Trans formal  inn  lm  [uc lined  Element 
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APPENDIX  A. 6 

STIFFNESS  MATRIX  (SYMMETRIC  AND  ASYMMETRIC  LOADING) 


Since  tlu»  stiffness  matrix  for  symmetrical  loading  is  obtained  as  a 
special  case  of  that  for  asymmetric  loading  we  discuss  here  the  Fourier 
harmonically-coupled  stiffness  matrix.  The  circumferential  variation  of 
the  displacements  at  any  meridional  station,  s,  may  be  expressed  as 


u(s,0) 


(o) 

u  (s) 


+  \  u  (s)  cos  j9  +  y  0  (8>  sin  J0 

4*?  j-1 


<°>  V  < V  .<*> 

v(8,0)  *  v  (s)  +  2^  v  (s)  sin  j0  +  v  (■)  cos  J0 


(1) 


w(s,0)  = 


4*1 

( <> )  V'  (  4)  (  4> 

w  (s)  +  y  w  (s)  cos  J0  +  y  w  (s)  sin  J0 

4*1 


4*1 


where  the  unbarred  and  barred  coefficients  are  defined  as  the  "A"  series 
and  "B"  series,  respectively.  Note  also  that  the  circumferential  dis¬ 
placement,  v,  varies  as  an  odd  function,  and,  therefore,  this  coefficient 
is  associated  with  sin  jO  (odd  function)  in  ,!A"  series.  The  "A”  series 
harmonic  number  j  takes  on  the  positive  sign  whereas  the  "B"  series  har¬ 
monic  number  j  takes  on  the  negative  sign. 

Similarly,  we  may  write  all  the  midsurface  strain  components  in  the 
form  of  Fourier  series  as  follows: 
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<°>  "  (})  *  (j) 

(e.'o  ■  <v.  +  2 ‘Vo  cos  je +2  <v0  ,ta  )6 

(<0  *  (  J  )  “  (  1 ) 

('.A  *  <*<>>«  +  >  <*  >  *ln  J0  +  >  («  J  cos  j0  (2a) 

80  O  80  O  £ 4  BQ  O  £  80  0 


(o)  A,  <J>  ”  _  <  0> 

<V0  '  (*g)o  +  2  ‘Vo  c0*  16  +  2  ‘‘o'.  8ln 


JO 


J*1 


Bending  Strains  are 

(0)  * 


X  -  X 

a 


(«)  “  (j)  “  (j) 

+  y  xg  cos  jo  +  y  sin  je 


V 


(O  ) 

*  (,) 

CO 

.()> 

x«e  + 

2x>e 

iml 

sin  J8  +  2 

X80 

<o) 

“  <>> 

00 

(  J  ) 

+ 

<D 

2 x* 

cos  jO  +  y 

\ 

(2b) 


Substituting  "A"  series  of  (1)  in  (16)  of  Appendix  A.l  and  setting 
these  results  equal  to  the  "A"  series  of  (2),  we  obtaln^conslderlng  only 
linear  terms ,  the  following  relationships: 

For  the  meridional  Strain, 


( i> 


<V.  -  f. 


2  %>co8 19  •  2  ” 


(  j  > 
(  s) 


cos  j0 


djo 

ds 


<  J> 

du  (i) 
ds  W  ds 


cos  j6 


Oa) 


j-i 
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From  (2a) 


®  ( n 

(e  )  *  V  (a  )  cos  j6 

so  £  80 


)  1 


From  Eqs.  (3a)  and  (3b) 


,  '•>  “>  m 

‘Vo  "t’"  da 


(3b) 


(3c) 


Similarly  for  the  shear  strain, 

“  <J)  \ 

v  sin  j0l 
iml 


<V.  -  r  Ct  * 


CD 

(2 


*  <  j> 

^  v  sin  j8  sln(? 


4-* 


+  -* 

be 


/  ”  U>  x 
|  X  u  cos  j0l 


(i)  ( n 

dv  v _ 

,  ds  r 


sin  tp 


( 4 )  “ 

•02 


sin  J8 


(4a) 


(  4 ) 

<e  >  *  >  <’  J  8ln  .10 

8Q  °  80  0 


<  se’o  ■  is 


4-1 

(  4)  .  (4) 

V 


(4) 


sin  <p  -  j  — 
r  r 


(4b) 


(4c) 


Omitting  the  first  two  steps  for  the  rest  of  the  strain  expressions, 

(  4  )  (4)  (  4  ) 

u 

(5) 

(4)  ...  (4) 


(J)  v  u  w 

(•0>o  " J  7  +  7  sln  w  +  7  008  9 


U)  £x'  '  +  u<J)^.j  +  ^  * 


8 


ds" 


ds  da  ds 


(6) 


V i ■V'l  ’  i  y  <P-r^ 


( J) 

[j  ^  -  ~  ain  9  w  +  coa  ®  ^ 

a0  r  tJ  3a  r  3a 

1  <J>  (,n 

-  ~  ain  9  coa  Q  v  +  j  u  J 

<4)  i  fia  (i)  j  <4}  ^<1}  (J)  »a  T 

L  -  -  -  w  +  -  cos  ®  v  -  f-  +  u  f*  aincp  1 
8  r  lr  r  3a  3a  J 


For  dome-ended  shells  of  revolution  (r  ■  0) ,  the  strain  components 
are  obtained  by  taking  the  limit  of  the  above  equations  aa  r  approaches 


zero. 


<<>  <i>  <<)  Sa'”  (1>  S? 

<Vo  ‘  1U  (Vo  '(Vo  ’to  •*  * 

rm0  irO 


( J ) 

(«  )  *  lim  (  «)  »  lim  3s  r  3a 

80  o  eg  o  1 


(4)  (J) 

v  ain  «p  -  Ju 


(J)  (J)  (J) 

3rSx  +t£x  .is  ,ln  ,  .  ,  5*  >a 

Um  45_5s__IiZ _ is  3  is  •  -J  g 

r-0  il  8* 

da 


Note  that,  at  r  ■  0,  we  have 


3r 

ain  ?)  '"“*  1 
3s 


3a  L* 


■  (i)  u)  (4) 

J  v  +  u  ain  +  w  coacp 


(O  *  lim  (a  )  -  lim 

V  t-0  6  °  t-0 

r**0 


<4>  (4) 


+lr  .‘n  *  +  * .«  *  +  «£(-*  15)] 
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■  j 


m 


& 

ds 


<  J) 


<  i) 

w 


(ID 


Note  that  the  first  derivative  leaves  the  expression  still  singular  for 


r  •  0  and,  therefore,  the  second  derivative  is  required: 


Similarly,  for  the  "B"  series  part,  corresponding  pole  strain  dis¬ 
placement  relations  are  obtained  by  replacing  unbarred  quantities  by 
barred  quantities  and  j  by  -'j. 

Now,  introducing  the  displacement  functions  of  Appendix  A. 4  into 
these  harmonically  coupled  strain  expressions,  it  is  a  simple  matter 
to  cbtain  W  in  (6)  of  Appendix  A. 5  or  in  the  expression, 

(  j)  0  )  (1  > 

X  “GW  Q  9  +  2  W  jQ  0 


(15a) 
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in  which  m  ■  cos  <P  and  n  ■  sin  V. 

For  axisymmetric  loading  all  terms  containing  J  in  W  are  zero  and 

Y  -  G  W  g.  ©  +  Z  ¥  Q  0  (15b) 

Comparing  this  with  (2.40)  in  Section  1  we  note  that 

Ah„6  -  £  W9  (16i) 

(16b> 

Note  that  in  the  nonlinear  membrane  strain  term  can  be  determined 

similarly. 

With  (15)  and  (16)  substituted  in  (2.40  a,b)  of  Section  1  and  all 
stress  components  of  oa^  expanded  in  Fourier  series  the  linear  elastic 


stiffness  matrix  of  the  form  (2.45)  of  Section  1  can  be  derived  explicitly, 


where 


rV2\  <i>  (j> 

■hi  I  £  W  T  a T  D  a T  W  £ 


rddde 


( j> 

T  £  bT  W  Q  rd6da 


1  0  0  0  0  0 
a  -  |0  10000 
0  0  1  0  0  0 


(17) 


b  - 


0  0  0  1  0  0 
0  0  0  0  1  0 
0  0  0  0  0  1 


l-v  V 

1  *eV 


*s 

0 

Ve* 


o  Et,v0. 
G  0 

0  E0 


Here  integration  along  the  meridional  coordinates  can  be  performed 
using  Simpson's  rule. 


WJ&  >,\ ' ■a5-«OT^W^T»»»,»3g*r^E9C35p??S,,^PT»^awra»'^-',5 «**<? ■«A*0‘S^«\/.*v. 


APPENDIX  A. 7 


EQUIVALENT  NODAL  LOAD  VECTORS 


We  consider  arbitrary  asymmetric  distributed  loads  of  Pi(s,6), 
Pa(s>®)»  a“d  Pb  (s ,9 )  acting  on  the  midsurface  of  the  discrete  element 
In  directions  s,  9  and  respectively.  In  matrix  form,  these  loads 
are  expressed  by  Fourier  series  of  the  form, 


<o>  V**  <j)  <  j)  V  (j)  Jj) 

U,(a,9)  -  £,(s)  +  >  C  p(s)  +A  £  P(») 

.  .  ~~  Jut 


where 


£(s,0)  -  [pi (8,9)  ,  pa (8,0)  ,  Pb(B,9)  y 


CO8J0  0  0 

0  sinJS  0 
0  0  cosJ0 


slnjU  0  0 

0  cosj©  0 

0  0  slnje. 


Deleting  the  last  row  of  matrix  jS  In  (4),  Appendix  A. 4,  we  write 


6  ■8/9>6  6 
3x1  3x8  8x8  8x1  3x8  8x1 


(2) 
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where  £  Is  the  3x8  normalised  interpolation  function..  The  equiva¬ 
lent  nodal  load  vector  PN  in  (2.38a)  la  derived  from  (2.20)  for 
harmonically  uncoupled  part, 


<  o) 


(0) 

p  \fi) rdSda 


or  in  matrix  form 


where 


P3  «.«)*♦  rds 


m  ZTt 

(o>  I  r  <o)  .  ( o )  ( o ) 

P  -  I  j  *T  p  w  rdeda  -[Pj  ,  P§  , 

(.)  *,0*v«»  ,rl  (.> 

Pi  -  2tt  I  Pi  (s;  Xjrds  -  *0  I 
•'o  ■'o 

(o)  r6  <o> 

Pi  -  2tt  I  Pa  (a)  Xarda 

•'o 

(o)  f  *  < 

3  ■  2tt  I  pa 


(oL 
•Pe  1T 


P3 


( o 

P4 


<0) 

(a)  *3rds 


(s)  X4rds 


)  P  (o) 

*  2tt  j  fa  (8 

"x> 

.3)  r  (*  <o>  i 

5  “  2tt  I  Pi**)  Xg  rda  '  j  Pa  (8)  Xsrdal 

•'o  Jo 

(o)  fA  <o) 

a  *  2n  I  pa(»: 

u 

f*  (0) 

-  2tt  j  PS  (a) 


(s)  Xp  rds 


(o> 

P-p 


( o) 


X7  rds 


0 

re  (•> 


P«  -*nl  Pj(b)  k,rdi 


where 


lit 


h  *  la  -  1-  |  .  la  ■  1  -  3(f)*  +  2(f)» 

»*•«}-  2(f)3  +  <f>*3  ,  V.  “  1«  "  f 

1,  ■  3(f)*  -  2(f)3  ,  la  -  t[-(f)3  +  (f)3] 


For  harmonically  coupled  "A"  series  part, 


<J>  Pf2n  (*>  (j>» 

L  J0  l'  E(8>  ~ 


rd0ds 


«  •  <J>  (o) 

Here  the  components  of  P  are  the  same  as  for  P  except  for  2tt  re- 

( o)  < ) ) 

placed  by  n  and  px  ,  etc.  by  pt  ,  etc.  The  components  for  "B"  series 
part  are  the  same  as  in  "A”  series  part  except  for  j  replaced  by  -j. 
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APPENDIX  A* 8 

CALCULATIONS  OF  STRAIN  AND  STRESS 


With  nodal  displacements  available  as  a  result  of  solving  equa¬ 
tions  of  motion  or  equilibrium  we  can  calculate  strains  either  by  the 
standard  finite  element  procedure  or  the  finite  difference  represen¬ 
tation  of  strain-displacement  relations. 


1.  STRAINS 

1.1  FINITE  ELEMENTS 

Using  the  notations  Introduced  In  Appendix  A. 4,  the  mid-surface 
normal  and  shear  strains  and  bending  strains  are 


_<j>  <j>  <i> 

V  -  W  Q  0 
6x1  6x8  8x8  §xl 


where 


2  •  l>.  ••e  «e  *• 


and  j  Is  the  harmonic  number.  Also,  the  engineering  strains  Y  may 
be  written  in  the  form  (see  (S)  Appendix  A. 4), 


(j)  (j) 

V-(G  +  2)WQ0  (2) 

3x1  3x6  3x6  6x8  8x8  8x1 


1.2  FINITE  DIFFERENCE  PRESENTATION 

From  harmonically  coupled  expressions  for  strains  as  given  In 
Appendix  A. 6  we  write  the  finite  difference  analogue  as  follows: 


•stJ  t 

■U  5 


<ll  tu'1’  t«  <I>  <|>  , 

e.  *r‘’5  ■•*■•>  '■ 

1  ,  (1)  <li  ‘|)» 

®se  "  7,  '•r*v«  "  v"  8irfp«  "  JU«  J 


lp  <j)  <i> 

7  [  jv,  +  u,  simp,  +  w, 

r* 


■  co«$. 


e 

where 


<j)  <i> 

X.  •  P»»  <3d> 

<  i> 

2  Cj>  8i.<P.w.  <i>  8ln<p,co»ip,  <}>  ( j ) 

[JW  -  J  - - - +  cos^v.x  -  - - - v,  +  J<P,u,  ] 

r«  ”•  r« 

(3e) 

1  r  i2  (  ..  COS5?*  a  _j_m  l 

:  Lr»i  +  j  —  *•"  V,  -  Pt8in«P,  J  (3f) 

r»  r»  r» 


<J>  a(,) 

®l»p*i  ■  ®i»p 


())  C|) 

>  P+-  i  +  ®i » p 


0(4>  . 
wa»w.j 

<*> 

®a»  p 

<*>  J*] 

®»>pt;  +  ®»»p 

L 

1 

v»  " 

2 

C|> 

<  J> 

<l>  <1) 

®3  »P+1  * 

%  » p 

Q 

®4*pU  "  ®4»p 

L 

1 

*9% 

L 

<  J) 

(4) 

®4  »  p+  X  +  ®4  »  p 

4 

Vi 

with  *  denoting  midsurface. 

Writing  (3)  in  matrix  form,  we  get 
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-<i)  <i  >  <i> 

-  B  8 
6x3  8x1 


<j)  <4> 
Y  -  (  G  +  Z  )  B  8 

3x1  3x6  3x6  6x8  8x1 


2.  STRESSES 


and 


T 

6x1 


The  Fourier  aerlea  representation  of  stresses  Is 
n  A 

<«>  V  <j)  (}>  > 

a  •  a  +/  C  a  +ZjC  a 

~  ~  £  -  -  ~ 


(i)  _Cj> 

where  C  and  <3  are  defined  In  Appendix  A. 7. 
(4)  <i) 

and  o  ■  E  V  ,  etc. 


3.  STRESS  RESULTANTS 


The 


lnplane.^st 


res  a 


resultants  are 


°tCdC 


ob9M 


oQCdf 


(4) 


(5) 


i 


(6) 


The  stress  resultants  in  the  {-direction  ere  given  by 


aw  sin<p  im  simp 

Qs "  aT*  +T“ Ms  +  Tir-  Me 

.  6Mfl  2s in*  6M 

%-7aT  +  — «.e+-aT 


(7a) 


Once  again  these  equations  are  expanded  into  Fourier  series  for 
summation.  It  is  also  convenient  to  write  (7a)  in  harmonically  cou¬ 
pled  finite  difference  form, 


(J>  <j> 

(,)  2siap,  (,)  M8e.^i-Ma0.* 

-i<M  V  +  — -  (Mafl).  + 


(7b) 


4.  STRESS  RESULTANTS  THROUGH  LAYERS  OF  COMPOSITES 

If  the  shell  thickness  consists  of  n  layers  of  composites  we 
write  (6)  in  matrix  form, 

(8a) 

(8b) 

where  k  denotes  the  layer. 
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APPENDIX  A. 9 

ELASTOPLASTIC  ANALYSIS  PROCEDURE 


1.  GENERAL 

The  yield  criteria  and  coordinate  transformation  of  plasticity 
matrix  together  with  elasticity  matrix  are  discussed  in  Appendix  A. 2 
and  A. 5,  respectively.  The  yield  criterion  to  be  satisfied  is  checked 
in  the  local  coordinates  of  the  fibers,  which  requires  transformation 
of  global  strains  and  stresses  to  those  of  the  local  coordinates  be¬ 
fore  the  state  of  yielding  can  be  determined.  In  what  follows  we 
discuss  procedures  involved  in  static  loading.  The  viscoelastoplaatic 
analysis  for  dynamic  loading  la  discussed  in  Section  2.8. 

2.  INCREMENTAL  LOADING  PROCEDURE 

As  demonstrated  in  Section  2,  the  Incremental  equations  of  equili¬ 
brium  is  of  the  form 

(  • )  { p ) 

(Kn*  +  Kn„  )d0"  -  dFN  (la) 

or 

<t> 

Kmh  d0M  «  dFN  (lb) 

( • )  (  p ) 

in  which  KNh  and  are  elastic  and  plastic  stiffness  matrices, 
respectively;  and  d0N  and  dFM  are  incremental  displacements  and 
applied  forces.  It  should  be  noted  that  in  the  above  equation  only 
the  static  elaatoplastlc  behavior  is  considered  for  simplicity. 
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The  given  load  la  applied  in  increments  as  shown  in  Figure  A. 9.1. 
For  convenience  one  dimensional  case  is  illustrated  here.  For  a  typi¬ 
cal  incremental  load  the  incremental  strain  to  be  reached  is  shown  in 
Figure  A. 9. 1(b).  Point  C  is  reached  through  several  steps,  each  step 
consisting  of  solving  (lb)  with  updated  plasticity  stiffness  matrix 
K**  which  is  calculated  from  plasticity  matrix  The  above  pro¬ 

cedure  is  referrred  to  os  the  tangent  stiffness  matrix. 

An  alternate  approach  is  to  rewrite  (lb)  in  the  form 


(•)  (p) 

K**  d0"  -  dFN  -  K*n  d0" 


(2a) 


If  the  load  increment  is  sufficiently  small  it  is  possible  to  write 
(2a)  for  a  load  increment  step  (s)  as 

< • )  ( p) 

Kaa  deM(s)  -  dFN(s)  -  K*N(8)dB"(s-l)  (2b) 


This  procedure,  called  initial  stiffness  method,  is  illustrated  in 
Figure  A. 9. 1(c).  In  this  case  each  iterative  cycle  is  controlled  from 
initially  calculated  elastic  stiffness  matrix.  The  product  of  plastic 
stiffness  matrix  and  the  Incremental  displacements  of  the  previous  step 
serve  as  a  pseudo  plastic  load.  Although  this  procedure  is  simpler 
than  the  tangent  stiffness  method,  it  requires  more  iterative  cycles 
before  convergence. 

The  transitional  and  unloading  elements  which  turn  to  plastic 
from  elastic  and  plastic  to  elastic,  respectively,  must  be  treated 
accordingly.  For  an  element  whose  equivalent  yield  stress  9  is  less 
than  the  current  yield  stress  d  at  the  end  of  elastic  analysis  or 
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dFa  “>  dot 

dFj  «>  dOi 


d@j  dS^ 


(a)  One  Dimensional  Representation  of  Stress-Strain  Curve 


(c)  Initial  Stiffness 


Figure  A.9.1:  Procedure  of  Incremental  Loadings 


1 
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previous  Increment,  it  Is  likely  to  expect  faulty  Incremental  equiva¬ 
lent  stress  ddi  larger  than  actual.  It  la  easy  to  show,  however,  that 
correct  for  such  transitional  element  is  given  by 


.(p)  . 

<toa-E(p)dVa  +  omaj{  -  a  (3) 


where 

_( p) 

dYt  -  (o  +  dox  -  CTmax>/<E  +  E(  p>)  (4) 


For  the  next  cycle  this  element  is  no  longer  transitional  but  elas- 

.(») 

toplastic.  In  this  case  dYi  above  is  modified  to 


.(pi 

dY. 


(p) 


[dVj  +  (5  +  doj 


a  >/E] 
max 


E 

BfE(,) 


(5) 


-<p> 

which  is  to  replace  dYj  in  (3). 

For  unloading  element  do  becomes  negative.  In  this  case  we  simply 
set  8*^*  equal  to  zero  and  only  the  elastic  properties  are  used  in 
the  next  cycle. 

Whether  tangent  or  initial  stiffness  method  is  used,  the  proce¬ 
dure  begins  by  scaling  displacement  and  stress  in  order  to  create 
impending  yield  (the  elastic  load  limit  of  the  structure  under  the 
given  loading  pattern).  After  the  elastic  limit  is  found  status  of 
yielding  is  checked  for  each  element  and  the  plastic  stiffness  matrix 
is  developed.  The  incremental  displacements  are  calculated  until 
convergence  with  updated  plastic  stiffness  matrix.  The  acceptable  per¬ 
cent  convergence  «  may  be  defined  as 


c  =  [{E(cr  +  djj)*  +  E( 5  + 


An  acceptable  convergence  la  conaldered  to  have  been  reached  if  c  la 
email,  aay,  e  5  57.  -  10%.  A  new  load  Increment  is  initiated  and  the 
above  steps  are  repeated  until  application  of  total  load  is  completed 


3.  INTEGRATION  THROUGH  THICKNESS 

( ») 

Calculation  of  plastic  stiffness  matrix  KNH  is  performed  by 

*(k  ) 

summing  the  Integrated  plastic  matrix  D  through  the  thickness  of  the 
shell.  It  is  possible  to  write 


Since 


"  K 

DC<*C  “  >  D  ( 


(h(tt)“  h(n-0 


K  )  »  2 

2.  (h<0  •  h(K-i)>/2 


J°c*dc  ■  ^  6<lt>(h*«:>  "  «*-»o/3 


we  obtain 


h(*)- 

(h?u )“  h(R.i)/2 


<h<«>-  hfm-t>/2 

W  a  rdOdf. 

3  3 

0»<k)  -  h(K.l)/3 


Here  h  is  the  thickness,  (k)  represents  the  layer  number,  and  §>  and  W 
arc  defined  in  Appendices  A. 4  and  A. 6,  respectively. 
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APPENDIX  A. 10 

DIRECT  NUMERICAL  TIME  INTEGRATION  SCHEMES 

1,  GENERAL 

If  the  dynamical  system  is  linear  and  has  a  simple  geometry,  the 
mode  superposition  method  may  be  used,  in  which  the  forced  response 
for  each  mode  is  calculated  by  way  of  the  Duhamel  integral  or  any  equiv- 
alent  method  and  the  total  response  is  obtained  through  superposition. 
This  approach  is  especially  attractive  if  low  frequency  bands  of  exci¬ 
tation  dominate  the  applied  loading.  Even  for  a  large  system,  conden¬ 
sation  or  component  modal  reduction  schemes  can  reduce  the  problem  to 
manageable  size  without  significant  loss  in  accuracy,  provided  the 
applied  loading  has  this  low  frequency  domination. 

When  high  frequency  excitation  is  significant,  however,  the  cou¬ 
pled  equations  of  motion  of  the  system  can  best  be  solved  by  direct 
step-by-step  integration.  However,  choice  of  inadequate  integration 
operators  result  in  unbounded  solution  or  unstable  solution.  In  the 
following,  we  discuss  various  integration  operators  and  the  solution 
stability. 

2.  FINITE  DIFFERENCE  OPERATORS  AND  ERRORS 

The  basic  differential  equations  of  motion  are  expressed  in  the 
form  of  recurrence  matrix  of  finite  differences  to  solve  dynamic  re¬ 
sponse  of  structures  based  on  the  following  assumptions: 

(l)  The  continuous  lapse  of  time  during  the  motion  is  divided  into 
a  series  of  small  finite  and  equal  time  intervals.  Within  each  time 


interval,  the  motion  of  the  structural  system  can  be  described  by  an 
ordinary  linear  differential  equation  of  motion. 

(2)  The  force  excitation,  or  the  displacement  excitation  applied 
to  any  part  of  the  structure  can  be  evaluated  numerically  at  any  de- 
signated  time. 

Numerical  Integration  techniques  fall  into  three  rather  broad 
categories  *  explicit,  implicit,  and  predictor-corrector  techniques. 

1.  Explicit  method 

a.  Forward  Euler  formulas 

b.  Runge-Kutta  formula* 

c.  Open-end  integration  formulas  (.predictor  formulas) 

2.  Implicit  Method 

a.  Backward  Euler  formulas 

b.  Closed-end  integration  formulas  (corrector  formulas) 

3.  Predictor-correction  method  (mixed  explicit -implicit  Method) 

3.  PARABOLIC  COORDINATE 

The  choice  of  the  finite-difference  equivalents  directly  governs 
the  accuracy  and  performance  of  the  procedure.  Consider  the  equation 
of  motion, 

M§  +  $5  -  F 

Introduce  a  small  time  interval  At, 

At  -  ta  -  ta_; 

Writing  §  for  a  parabolic  variation  of  ®  within  the  time  Increment, 
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•  •  1  _ 

0  m  I  -  I  r®  L  m  -  2®  +  @  .  \ 

xt  (At)aL~"fl  xn-1  > 

Substituting  in  the  equation  of  motion 

H  ■(It)5’  C&»i  -  2®n  +  0  ..tJ  +  Kg.  -  £ 

®»+l  -  M"lF(At)*  +  [2  -  M‘lK(At)»3  ®#  -  ®  #.l 


(1) 

(2) 

(3) 


Figure  A. 10.1  Parabolic  Variation  of  Displacement 


For  free  vibration  with  single  degree  of  freedom 

®«n  *  C2  ’  jj”(At)*]®B  -  0^  (4) 

Introducing 

0,  "  AP"  (5) 


where  A  is  the  arbitrary  constant  to  be  determined  from  initial  con¬ 
ditions  and  P  is  a  number  to  be  so  chosen  that  (4)  is  satisfied. 
Substituting, 


Ap«fl  +  C~<At)«  -  2]APB  +  Ap“-1 


0 


Dividing  through  by  AP""1 

P*  +  C^(At)»  -  2]p  +  1  -  0  (6) 

Consider  various  ranges  of  At 

(1)  0  <  At  <  2/- 

First  investigate  At  -  ^2/^  «  .225(2tt)/^  so  that  At  -  .225  T 
2t? 

where  T  -  jjjr  ,  natural  period  of  the  system.  In  this  case  Eq.  (6) 
becomes 

P*  +  1  -  0 

P  *  +•  /-I  *  -  i  *  e  +  y~ 


(7) 

(8) 
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to  ^  . 


f' 


As  At  changes  from  0  to  2  ,  the  factor  In  (10)  varies  from  1 


(2)  At  > 


Consider  first  At  *  3 


7f  •  e«- 


(6)  becomes 


P*  +  7?  +  1  -  0 


(ID 


giving 


-  -.1459  -  -e 


-1.925 


P  *  -6.8541  ■  -e 


1.925 


Substituting  these  into  (5), 


0n  -  A(-l)”  e'1,925w  +  B(-l)n  e1,925" 


A  'cos  n  tt  sinh  i.925n  +  B/  cos  n  tt  cosh  1.925n  (12) 


Substituting  n  -  —  ■  t  J  m 


0.  ■  [A'  sinh  . 642  t  /-  +  B '  cosh  . 642  t  I— 
*  •*  m  J  m 


]  cos  1.047  t  /— 
1  m 


(13) 


The  primary  effect  of  numerical  integration  methods  is  to  Introduce 
hyperbolic  functions  of  time  in  the  arawer.  Since  these  functions 
Increase  Indefinitely  with  time,  the  result  is  divergent. 


For  all  values  of  the  time  increment  At  ;  ,  divergent  solutions 


of  the  time  (13)  result 
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4.  CUBIC  COORDINATE 

Writing  the  acceleration  and  velocity  for  a  cubic  variation  of 
0  within  the  time  increment, 


Figure  A. 10.2  Cubic  Variation  of  Displacement 
Substituting  into  equation  of  motion 

M  (a?)*  [2£»  ’  5®«-l  +  +  l£»  *  l  (15) 

or 

(2  +  M"lK(At)»)0a  -  M"lFAt*  59^.1-41^.,  +  tf,_a  (16) 

For  free  vibration  with  a  single  degree  of  freedom 

[2  +  “(At)a]0#-  58,.;  +  48,.*  -  9 ,.3  -  0  (17) 


Substituting  (5) 
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[2  +  -  A t* J Ap  n  -  5Ap»"1  +  4APa"*  -  Apn’3  -  0  (18) 

tn 

Dividing  through  by  A8”"a 

[2  +  -(Atf]^  -  5^  +  48  -  1  -  0  (19) 

m 


In  this  case,  the  errors  are  a  small  drop  in  natural  frequency,  ,961 
in  place  of  1  due  to  the  presence  of  decaying  exponentials.  In  the 
limit  when  At-O,  the  coefficients  2,32-*  »  ,  ,0095  -♦  0,  .961  -•  1.  The 
solution  approaches  the  exact  solution  as  At  approaches  zero. 


Case  (3).  To  judge  the  errors  for  large  values  of  At  consider 


At  -  98 J 


■11805, 


•  -  ^umf,  > 

x(B '  sin  .I83t  ]&  +  c'  cos  ,l83t  I®  ) 
1m  f  m 


The  errors  are  a  large  drop  in  natural  frequency,  .183  in  place  of  1. 
As  At-  *  we  have  coefficients  .805,  .1424,  and  .183  -  0.  It  is  seen 
that  large  errors  may  result  for  large  values  of  At,  but  they  will  be 
always  convergent. 

5.  GENERAL  FORMULAS 

Consider  a  time  interval  At  and  assume  variations  of  acceleration 
to  be  constant,  linear,  of  a  step  function,  or  of  any  other  form.  For 
example  let  us  take  a  linear  variation. 


128 


«  -  Si  +  i)  fe 


By  integration, 


&  *  Sat  +  (8*x-  |n>2^+  *1 


where  Ci  *  0„.  By  substituting  t  *  0 

»•§»  t+&*»  -®«)24t  + 


By  further  integration, 

®  "  ®n  +  "  ®n>  6AT  +  St  +  c» 


where  c»  ■  0_.  By  substituting  t  ■  0 


..  ta 


(20) 


®  “  %~2  +  §«>6Al  +  ®»t  +  ®» 


At  t  ■  At  we  have  <§  -  0nf  1 ,  ®  ■  ®  n+  i 


&*l  "  ®*At  +  (^n  -  ®.)  f*  +  §  -  6,  +  ^  +  g* 

2  .  -  &rf  +  &m  -  §,)  ^  +  i  «  +  9, 


*  2.4C  +  5  2.  (At)*  +  I  g»,(At)» 


Newmark  suggested  the  following  general  form  of  and  0^;  by  intro* 
duclng  parameters  Y  and  8  which  characterize,  respectively,  artificial 
damping  and  patterns  of  acceleration  between  the  time  interval: 

t,n  *  L  +  (1  *  Y>  +  ig*xAt  (21) 


(21) 
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®»U  -  h  +  g/t  +  (|  -  8)0,, (At)*  +  P0„.i(At)*  (22) 

For  Y  ■  0  produces  negative  damping 

Y  >  ^  produces  positive  damping 

Y  =*  £  §U+i  *  g»  +  \  At(©n+t) 


Figure  A. 10.2  Step,  Constant,  and  Linear 
Variation  of  Displacement 

It  can  easily  be  seen  that  we  must  choose  Y  ■  ~  to  avoid  artificial 
damping  and  8  »  £,  8  ■  and  8  *  ^  correspond  to  accelerations  of  con¬ 
stant,  linear,  and  step  function,  respectively. 

It  is  interesting  to  note  that  Eq.  (21)  resembles  the  truncated 
Taylor  series  expansion  with  particular  choice  of  coefficients, 


®*+l  "  k  + 


(23) 
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Substituting  these  into  the  equation  of  motion  for  the  time  step  cor¬ 
responding  to  nfl, 

1  +  +  ®.At  +  &  +  -f  -  g^At*)  -  Fj,+  i  (24) 

or 

<a+!^r)|»i  -e*»  -  te.  +  a^t  + 1.  (25) 

From  this  we  calculate  and  0^,  subsequently. 

6.  RATE-DEPENDENT  PROBLEM 

If  viscosity  is  present  the  general  form  of  equation  of  motion 
becomes 

MG3+C0  +  K0-F  (26) 

where  C  is  the  viscosity  matrix. 

Considering  a  constant  acceleration  and  substituting  (22)  and 
(23)  into  (26)  for  the  nfl  step, 

M0.H  +  C 

+  K  (0.  +  At  %  +  f£-  0nn  )  -  P,a 

<l£+^£  +  4r£>®»n  -  Em  -  Sm  <27> 

Em  -  £<®»  +  4  ®»)  +  k(®«  +  Ate,  +  ^  ®») 


where 
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Recursive  equations  are  then  (27) ,  (23)  and  (22)  to  be  solved  repeat¬ 
edly  for  all  time  increments. 

7.  QUASI-STATIC  PROBLEM 

In  the  absence  of  inertia  force  we  have  the  equation  of  the  form, 

C  0  +  K  0  -  F  (28) 

Either  the  displacement  or  displacement  rate  may  be  assumed  to  vary 
linearly  within  a  small  time  increment.  Both  cases  are  considered 
below. 

7.1  LINEAR  VARIATION  OF  DISPLACEMENT  RATE 
For  a  given  time  Interval, 

&•  &  +  s.) 

S*  !»  +  24t  ®*f»  "  ®*>  +  ®" 

For  t*At 

+  r<V  •  £*>  +  9- 

or 

9»+i  +  (29) 

Substituting  (29)  into  (28)  for  the  «+i  step 

£  tu  +  i  +  ¥  K 

or 
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<£+r  -£*> i-&»i  <30) 

where 

-  k<».  +  Is-  s»>  <31> 

Here  (30)  and  (29)  constitute  recursive  equations. 


7.2  LINEAR  VARIATION  OF  DISPLACEMENT 


We  consider  a  linear  variation  of  displacement  so  that  the  aver¬ 
age  displacement  at  the  mid-interval  is 


-  (6..1  +  8,)/ 2  (32) 

The  displacement  rate  at  the  mid- intervals  is  approximated  as 

®..4  '  <8„i  -  e„)Mt  (33) 


In  view  of  (32)  and  (33)  and  rewriting  (28 )  for  the  mid-interval  , 
we  obtain 


28. «-  20. 


c(- 


At 


■)  +  K  ®n+..  -  Zt»\ 


or 

<ft£+S2^j*».*jk  +  ft  £8. 


using  (32)  once  again, 


-2<ft  £+«<E»s  +ff£g»>  -s.  (34) 

which  may  be  solved  repeatedly  to  any  extent  of  time  desired. 
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APPENDIX  A. 11 
MASS  MATRICES 

1.  GENERAL 

There  are  two  types  of  mass  matrices  In  use:  lumped  mass  matrix 
and  consistent  mass  matrix.  The  lumped  mass  matrix  is  a  diagonal 
matrix  contributed  by  the  tributary  area  around  the  node  with  equiv¬ 
alent  mass  (weight  divided  by  gravitational  acceleration).  The 
consistent  mass  matrix  is  derived  in  Section  2  and  given  by  (2,39) 
in  terms  of  the  normalized  interpolation  function.  Although  the 
lumped  mass  matrix  is  simpler  we  use  the  consistent  mass  matrix  in 
the  present  study  because  computer  coding  does  not  present  any  special 
difficulty. 


2.  THIN  SHELL 


The  general  form  of  consistent  mass  matrix  (2«39)  is 


M  /•  2rr 

-  I  I  Pt 

Jq  Jo 


1ttHr(s)d0ds 

N 


or  in  matrix  form 


M  ■  2ttp J* 


*>)  £(b)  rUQds 


or 


M  *  2ttp  2,t  t  2. 
8x8  8x8  8x8  8x8 


in  which 

r* 

L  "I  ST  S  r<s)ds 
8x8  8x4  4x8 


(1) 


(2) 
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4 


Here  and  are  given  In  (5),  Appendix  A. 4  and  the  Integral  L 
la  easily  obtained  by  computer  program  using  Simpson's  rule. 


3.  THICK  SHELL 

The  consistent  mass  matrix  corresponding  to  a  thick,  shell 
without  bending  represented  by  the  plane  strain  element  described 
in  Appendix  A. 12  can  be  obtained  similarly  as  in  the  previous  sec¬ 
tion.  We  have 


ing  by  means  of  Gaussian  quadrature  is  quite  simple 
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APPENDIX  A. 12 

THICK  SHELL  ISOPARAMETRIC  ELEMENT 

1.  INTERPOLATION  FUNCTION 

Consider  the  four-sided  element  with  nodes  1,  2,  3,  and  4  charac¬ 
terized  by  two  coordinate  systems  r-z  (cartesian  coordinate)  and  ?-l] 
(isoparametric  coordinate)  as  shown  in  Figure  A. 12.1.  The  isoparametric 
nodal  values  are  either  1  or  -1  measured  from  the  origin  of  §-i|  coordi¬ 
nate.  It  is  possible  to  write 


r*  5.1  (1) 

where 

r  -  [  i  ?  'i  sn  ] 

a  -  [  al  a?  a3  a4  3 


Here  a^,  etc.,  are  constants  to  be  determined  writing  the  nodal  values 
of  X  by  substituting  the  nodal  values  of  5  and  T1  we  obtain 


r  ■  A  £ 

where  _r  •  [  ri  r,  r3  r4  ]T 

"  1  -1  -1  l 

1  1  -1  -1 

A  * 

~  1111 
1-1  1-1 


(2) 


from  which 
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Substituting  the  above  into  the  expression  for  r  gives 


r  -  C  l  5  »l  SH  J 


1111 


-l  l  l-l 


-l-l  l  l 


l  -l  l  -i. 


where 


r  "  1  £ 


ti  -  !<i-S)<i-l]> 
♦a  -  |(1+?)(1'1|) 
♦a  »  £<1+?)<1+H) 

*4  -  J<1-5)(1+T1) 


Similarly  z  can  he  written 


Here  ^  is  called  the  isoparametric  interpolation  function. 
Combining  (3)  and  (4)  we  can  also  write 


In  matrix  form 


2  -(l]*!2 


where 


Xt  -  tlrX' 


i- ill’ll 


♦a  til 

ta  tsj 


ti  0 


0  ta  0 


ti  0  ts  0  ts 


X-C  rx  r„  ra  r4  zx  z9  za  z4  ] 


X  -  [  rx  zx  ra  z8  Xq  z,  r4  r*  ] 
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* 


4 


Here  1,  2,  3,  and  4  Indicate  the  node  numbers. 

It  Is  now  assumed  that  the  displacements  u  and  v  are  also 
related  by  the  same  interpolation  functions  such  that 


where 


u  *  f  u,  V  ■  ♦  V 


or 


«  "  Cui  us  uJt*  l  m  C vi  vs  va  v43t 

u  -  J^d-*  -7i  +-n))uj  +  (i  +  * 

+  (1  +  t  +  T)  +  fT))Ug  +  (W  +  T|-  ^l)u4l 
-  i  I  d  +  Bn  ?  +  Cn  T)  +  *l))uB 


Mail 


B  -  -1 

l 

B  «  1 


B  -  1 

i 

B  -  -1 

4 


C  -  -1 
1 

C  _ 


C  -  1 

3 

c  *  1 

4 


1 

D  -  -1 
a 

D  ■  1 

3 

D  -  -1 

4 


Similarly 


8 


v  ■  i  j]  (1  +  B  »  +  C.  1]  +  D.  rT|)V/ 

B»|  *  '  "  "  * 


(5a) 


(5b) 


Since  r  is  given  by 

*  -  i  C(l-u)<l*Tt  -*x+  (l+f)(l-T))rf  +  (1+^(1^))^  +  (W)(U-r|)r4] 

or  t  -  %  [,  R  T  +  (R*)?+  (R'g)T^  (R'C)<J]  ]  (6) 


Rt  “  (*!  +  *m  *  r3  +  V 

«.  •  <->,  +  ».  +  v  v 

K.  *  <-rt  *  C.  +  r.  +  V 
«c  *  <V  +  VC<) 


where 
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Similarly  t 

*  -  *<Qt  +  Q„1)  +  Q^)  <7> 

Where 

0,  -(*,+&,+  +  a4) 

Qa  -  (-*!  +  a„  +  z3  -t4) 

Q„  -  (-k4  "  »,  +  *,  +  *4) 

%  *  <ri  *  r,  +  r3  -  r4) 

Substitution  of  interpolation  functions  into  the  strain-displacement 
relations  will  lead  to  partial  differentiation  of  these  functions.  Thus, 


or 


This  gives 


in  which 


&• 

AI  +  A*- 

.  J5S 

hr 

ft? 

5s5 

5? 

hi  . 

m 

X£  4.  hi 

21* 

Wl 

hr  * 

Wl 

5* 

*  61) 

A£ 

Aif 

ht" 

he 

K 

he 

6? 

hr 

6? 

At 

AT 

Al 

i* 

-  i 

h* 

61) 

»  &7j 

6T) 

_5*_ 

_6z_ 

it 

'A*  _  'AS" 

hr 

6H  6? 

6e 

"  A 

A 

-M  2lX 

At 

_6*_ 

61)  6* 

■» 

1 

hi) 

A  ■ 

•  A£  _  A 5 

ail 

u 

U' 

6T)  h? 

WlJ 

Here  J  is  called  Jacobian  matrix, 
tion  required  in  (8), 


(8) 


(9) 


(10) 


Performing  the  partial  differentia- 
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ff  "  *[  <-r,  +  ra  +  r3  -  r4)  +  (r  -  r,  +  r3  -  r4)i|] 

•  Similarly 4 

-  if('^  +  z„  +  ^  -  z4)  +  (Z,  -  zf  +  Z3  -  Z4)f)] 

^  “  4  C("*i'ra+r3+r4)  +  (rj,  -ra+r3  -r*)5] 

•jr  *  J  C(“zi"z3+z3+Z4)  +  (Zl"2«+*3-z*)§3 

?iM  «L  [(*r1+r,+r;1-r4)(-zl-z.+z3+z4)  + 

t>?  tf|  16 

<r1"r»+r3"r4><"Zl"ZB+Z3+Z4)l)  + 
(-r1+ra+r3-r4)(zl-zB+z3-Z4)(r  + 

<ri  "re+r3  "r4  ^zi  ”za+z3  ’z4  >  0l3 

And 

"  16  C  <-ri’r.+r3+r4><-Zl+Z.+Z3'Z4)  + 

<’ri  ‘r*+t3+r4  >  <zi  "Z«+Z3  “Z4  >1)+ 

<ri  ■r«+r3  'r4  >  <Zl  ~ZS+Z3  ”Z4  >  ?T}3 

Substituting  these  into  Equation  (8)  and  denoting  that  r(,  -  r,  - 

A  ■  8  C  (Z4  SZ13  ”Z13  Z43  3  +  ^r34*la"ri3Z34^  +  ^r4  1  Z*3  "r»3  Z4 1 

Similarly, 

“  ^(“Vj +*.+*:, -v4)  +  (vt-vp+v3-v4)T)] 

■  ^[(-V^V  +^+V4  )+  (v1-V3fV3*V4)|  ] 


(11®) 

(Ub) 

(lie) 

(lid) 


etc., 

>1  3 
(12) 

U3«) 

(13b) 
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M  M  -  ll  ^  "v‘  +v*+v3  "v*  X "  ri  *  r9+ri+r* )  + 
(-Vj+v.+Va  -v4)  (tj  -r,+r3  -r4  )  $  + 

(Vj  -  va4va  -v4  )(-r1  -  r#+r3+r4  )T)  + 

(vt  -v»+v3  -v4 )  (rt  -r*+r3  -r4  >51} 


Hl.hL  m  m  M 
M|  ft?  6?  M| 


*  8  C<vi*«*iri3^«*M*4%i>  + 

(Vir34+V3r43+V3*3l4V4ri*>?  + 
(Vir83+V31P4.!'4V3rf4'H4*3»  >1)  3  ' 


(14a) 


^  ^  "  16  {(■v1+v,^3-v4)(*  -zB-ft3+*4)  + 

(-Vj+Vn+Va^Xz^Bp+tg-^)?  + 

(v,  -v,  -V4  )  (-*,  +*4  +*4  )T|  +  (vi  -vf+v3 -v4  )  (rt  <-3s  •*%  )Sfl3 

(vl-vt+v3-v4)(-«,+*g+*9-*4)^  + 

(-vt  -v#«hr3+v4)(*j -z#+z3  -*4)’I)+(vi  -v*-h&  -v4)  (tj.  -e*+*s  -*4  >5^3 

Ch«r«£or« 


ZH  M  .isil  -  i  f(v.  *94+Va*3|  +va,!4fl+v4*  13  3  + 

6S  6T|  Ml  6'  *  l'  1  >4  9  31  3  43  4  13 

(V,*43+Vf*34+V3*|«+V4*JM)?  + 

<V,  S«+Vl4+V3*4»+V4**l>Tl}  <14^> 


We  will  make  use  of  equations  ( 14a  >  and  ( 14b  )  in  evaluating  ^ 


In  terma  of  Ti»vg,vs,v4  and  nodal  coordinate*  (r,,*,),  (r3 ,«3) 

and(r4,*4). 

From  equation  (.  ’)  '■)  we  can  write, 


hi  .1  +  h*  .  Ml  (15a) 

a*  4  L  wi  &s  as  wU 
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and 


3SL 

.If  bz  JSL  .  a*  in] 

A  L«r)  *  h*  5r  571J 

(15b) 

Using  equation  (  14a.)  and  (14b  )$ 

m  . 

6* 

8A  fvl  r4*+verla+v3 ^4^4 ra i  >  + 

<V1  r3  *+V.r4S+#r3  r»i+V4  ri  *)  ?  + 

<vi  r«,+VBlf4|+Va '14^4  r3a>*n} 

(16a) 

and 

^  *  ll  f(  1^3*48^4  Z13>  + 

<Vl  «48'tVa*S44V3  ZJ  s+v4za  1 )  # 

<V1  Z3  a+v?zl  4^3  Z4  1+V4  za3  'T|} 

(16b) 

Denote 

0*  k 

M! 

-  n(r4s+r34f»-ra3Tl) 

M« 

"  n(rj3+r43?+r4^> 

H, 

*  n(rP4+r«irfr!4T]  ) 

M4 

"  n(r3l+rl??+r3aT)) 

h 

-  n(sa4+z43^z3PTl) 

h 

-  n(*31+*34*tfzj4T)> 

• 

h 

*  n(z48+zu^z4^) 

h 

■  <X*ia+*#i  f^wTl) 

Substituting  these  into  equations  (16a)  and  (16b)  yields 

O' 

6* 

-  hr  v 

(17a) 

OL 

•  L  v 

(17b) 

br 
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where 

M  -  tM,  M,  H,  MJ 
L-U,  I,  l-,  Lj 
and 

v  <Vl  va  v3  v4  ] 

Equations  (i-17»a)  and  (1-1.7 b)  can  also  be  written  as 


fr  ■  $n  <*„«,. #c, ,ii)», 

'  Cal 


where 

Ar, 

at 

r4’r3 

Bri 

- 

r3 

Cr» 

■ 

Vr3 

Ars 

- 

ri“r3 

Br, 

m 

r4  ”  r;< 

<*• 

■1 

V*l 

Ar3 

m 

rp"r4 

m 

r«-r! 

<*3 

m 

ri’r4 

Ar4 

w 

r3“rt 

Br4 

"rrra 

Cr4 

m 

r3-^8 

and 

Azt 

m 

23-24 

BZj 

m 

Z4*Z3 

CZj 

m 

Z3*Z^ 

Az» 

m 

Z3-Zl 

Bz. 

at 

Z3  ”Z4 

Czs 

m 

*1  ”Z4 

AZ3 

ar 

Z4"ZB 

Bz3 

ai 

Z1  “za 

CZg 

m 

Z4"Z» 

AZ4 

m 

Zl'Z3 

BZ4 

- 

Z2_Zt 

Cz4 

m 

Z8*Z3 

Similarly, 

0  <Wl 

“  "Ain<4'"+B”' 4C"11  )u" 

Using  equations  (  18a, l),  (19»,b) , (5a .  .)  and  (  5b  )  the  following 


(18a) 

(18b) 


(19a) 

(19b) 


relations  can  be  obtained 


144 


fiv  &v 
ftz  '  &Z 


S  *  fff  Ar.Ar»+(Ar.Brn+Br.Arn)^pCr»+Cr.Arn^ 
(■  1  »»1 


<Br.Cr«+Cr.BrnM  +  +  +  C„Cr  ^Jv*. 

5  f  -j  i: na  t4”A”+  <*»*..  +  +  <A,.C,.  +.«,A.>’i + 

S  S  ■  J  £  tf  f  *..*.«  +  +  +  c, .*..)!  + 

+  *,.C,.)PI  +  +  C,.C,.T|^1»  .V, 

%  f  '  !  !  »  {A,  A  ,+  <»,  „+  ",  A  .> •*  <*,  .8,  +*,  ,8,  ,>D  + 

Or  &Z  *  l  if*  1 

®,.c,.+  e,,*,,)^  +  b,.b„^  +  e„c„T|"l*.u, 

j£  J'nA  *-.+  <8,.+A,.®„  )'  +  <8,.+A,.8.  )1|  + 

(C  B  +B_  C  +A  D  )tfl  +  B  B  f"  +  C  C  n»  +  B  D  *»ti+ 

v  r«  n  r i  n  r  »  n  '  "'I  r»  n  '  r  »  n  'I  r  ■  n  r-  'I 

Cr.Dn  P]*}v.un 

57  P  -  \  s  &(A,.+  <A,A  +B,.)?+  «,.C,  «,,)!,  + 

■■1  nwl 

«..c.  +A,A  )w  +  >,.8.  f  *  8ItC,  D  +  B„D,  fTr*C„D.  «f)Y,U, 

?§[  ‘  *.  fe  f*”+  (,-»+,»A-.)«+  (C.  \ .+  °r. Ml  + 

**  1*1  0*  * 

<».  8„+C„  B,+D.  \„M  +  +  C.  C„  Tf  +  D.  B„.?T|  + 

D.  c,  ."I* I  u.u, 

r  S  '  *  t  ^  tA,  <8.  A.,+C,.)T|  +  (B.C..+ 

8.  +  »«»,.?  +  8.  C„t|  +  D.  B,„,  T|  +  D„  C,.<n  )«.«, 
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4  4 

“7  “  S  Z  Cl  +  <B*4Bb  )?  4  (Cm+CB  s)tl  4  (Dm4D,  4  B*C,  4C*  B,  )gT|  4 

r  r  ■*!  a«i 

B„  B,  £c„  CB  TH-  <»,*,  4Dm  Bt  >W  4  <CB  D,  4  DmC,  )5^DB  D,  (i\  } 

2.  STRAIN  ENERGY  AND  ELEMENT  STIFFNESS  MATRIX 
Strains  In  the  element  are  given  by 


- 

"bv 

Yz 

bz 

bu 

Yr 

m 

br 

Ye 

u/r 

v 

bu  bv 

Y 

zr 

bz  +  br 

-  „ 

Element  stress  strain  relationship  can  be  written  as 


! - 

N 

®11  °12  °13  °  " 

1 - 

N 

I5" 

U 

r 

°21  °22  D23  ° 

\ 

■r 

ae 

°31  D32  °33  ° 

a 

GOOD., 

y 

rz 

44 

rz 

■»  m 

M  <**> 

i..  which 

Djl  "  CEj  (E|.  _Et^tl  ) 

D12  »  CET(ELvTT4ETv*r) 

°13  -  I>23  "  CEtElvtl(1+vtt) 

^22  *  CET(Et,-ETvTu) 

D33  -  CE?(i-v*T) 

°44  "  Gt 

C  "  [Et(1-vtt)  -  2Et vTt  (14vtt)] 


aaaairitta  as. 


zrCTzr)dv 


Strain  energy  in  the  element  la  given  by 

u  •  I  A'2.,av 

•?///(V*+vro*  +  V,+  v; 

■  iff f  VDll*z  +  »12\  +  D13vo  +  VD2lVD22  VD23V 
2  +  V^lVD32VD33VVVr*>l  dd'ldr<"! 

l*"*'  //drdz=jjT|j|d5dT| 

u  .  I  /2rr  ^  f  f/D  &  to  +  D  M  ,  n  + 

2  y  J  *•  n  6z  °12  ftz  &r  °13  ftz  r 

(D  iEM+n  i&hii+n  M  u\  + 

*1  hr  P  +  °22  *r  23  $F  F  + 

‘  u  ftv  ,  u  ftu  u  u  v  . 

°31  r  hz  +  °32  r  *r  +  °33r  r  ^ 

,n  f  &  0°  X  ftu  ftv  J.  &•  &U  n  r  I  T I  Aniitr 

^D4A^  hr  hr  +  r>r  ftz  +  ftz  ftr  +  &z  &z  ^  lJldT1d^  (20) 


•1  -1 


1  1 


?.n  mjj  f  fAriAr.+  (Ar.B,B+»r.Arn)^<Ar>Crn+Cr.Ar.)Tl  <Br.Cr.+Cr.Br.>CT 


-1  =1 


+  BraBr«V+  C,  .C,./)  rd.ld? 


1  1 


■1  -1 


(Br.C>.+Cr.B,^;'l  +  Br.BlK^+  ^.C^T]3}  r  dr,d? 


1  1 


■■  ’//  S  (  A. A „+<A, +  (A,.c,.+c..A,.>TH- 


■1  -1 


(C,,B,.+B*nCt.)^  B„BIn^  +  Cs„Ctn7)a  }r  d^d? 


1  1 


■  n 


■//  a  tA. A „+  AA.+“,.\, + 


-1  -1 


<»..Cr  .>01  +  B,.Br,?'+  *1M? 
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P-M,  “//fefA r«+<Br.+Ar.B«  >S  +<Cr.+Ar.Cn  )n  +  <Cr.B«+Br.C«  +A-.D»  >  SH 


-1  -1 


+  Brl«,  r*  +  C„C,  #  +  B„D„  *1|  +  C„D,  pf )  Mi; 


X  X 

Q„  -/Jfe  (*,.+  <A,.»,  +b„)f+  <A,.C„  +Ci„>T|  +(B,,C,  +C„B„  +A..D,  )P1  + 


-1  -1 


B*A  ?+  C..C„  Tf+  B,.Dn  ^  +  C,.Dn  5  f  W? 


X  x 

‘JJ\ 2  (*-«  +  A-.)F.+(C.  Ar.+C..>1l  +  <B»  °-.+C«  B-.+D»  *..>CT  + 


-1  -1 


Bm  Br„?*+Cm  Cr.jf  +Dm  Brn*^  Bm  Cf  n$Tf  }dT|dP,  =  Pn, 


S-»"/ /  32  fA*»+(BB»+Bn.  A,.)5+(CwAtB+ClB)n  +  (B-,Cfa+CBBt)1+Dm  Ar„)5n 


-i  -i 


+  »„  B,„?  +  C.  C1„T|=+  D,  B„  pi)  +  D„  C,,P|»)<1T|B?  »  Q, 


T.„  ’//nr?  f1+  <B«  +B"  )?+(C«+C.  >B+<D.-+D.  +B»  C.  «.B.  )HT1  + 
-1-1 

a 

Bm,B^a  +  C„C„ -V  +  <BmD»+Dm  B«  >5  D 
(CmDn  +DmC«  +  D.  D«  lJ!dl*dP’ 

These  integrals  to  appesr  in  (20)  -are  evaluated  by  using  sin  point 
Gaussian  quadrature.  With  the  notations  above  we  can  rewrite  (18)  in 
the  form, 


»  -  J  <2n>  E  E,  fBnF.,V.»,+D  2C„V.  u,+H,g>v.U.+D210,.u.*„ 

+  D22^nU«U"+D2^.nu.u»+D3lPrt»u.VD320».u.u. 
+  D33T.nu.u,+DA4(I..V.V«‘rt*«U'«V»+  -  •U»)) 

,44 

U  “  o'  (2n)  E  E  f\,v,v,+B,.v,u|«l,u|v+8I1,u|ui] 

4  i-tn-l 
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*•» j  i”. 


The  stiffness  matrix  K  Is  Identified  as 


4  4 
E"  5 


where 


;  fA.n  B./I 

1  1C.,  B.J 


A  -  D..F.+D,  ,1.  „ 

■  n  XI  * n  44  ■ n 

*  d12G..+d13P.»+D4*I‘.. 

C.»  ■  D2lG..+D31P..+VK. 

E»  ’  D22I..+D230..+D329-+D33T..+D44f.. 

The  stiffness  matrix  obtained  from  Eq.  (21)  Is  of  the  form, 


V1 

v2 

v3 

*4 

U1 

u2 

U3 

u4 

A11 

A12 

A13 

A14 

I 

1 

B11 

B12 

Bn 

B14 

A2l 

A22 

A23 

A24 

1 

1 

B21 

B22 

B23 

B24 

An 

A32 

A33 

A34 

1 

I 

• 

B31 

B32 

®33 

B34 

,  ^ 
i-1 

• 

• 

A42 

A43 

A44 

1 

l 

1 

B41 

B42 

B43 

B44 

cu 

C12 

C13 

C14 

1 

1 

1 

E11 

E12 

E13 

E14 

C2l 

C22 

C23 

C24 

1 

1 

| 

E21 

E22 

E23 

E24 

C3l 

C32 

C33 

C34 

1 

1 

E3l 

E32 

E33 

E34 

C4l 

C42 

C43 

C44 

1 

1 

E41 

E42 

E43 

E44 

Elements  of  stiffness  matrix  thus  obtained  .•  should  be  rearranged 
so  that  they  .  conform  to  the  applied  force  vector  (or  computed  dis¬ 
placement  vector). 
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3.  FORCE  VECTOR  DUE  TO  BODY  FORCES 

The  force  vector  due 


0  O' 

*3  0  *4 


The  components  of  force  vector  become 

F-[0  Flt  0  Fjx  0  F3f  0  F4i]t 


in  which 

Fw" 

-2npt  J1  rl(l-?-il+?n)r|j|d1]d? 

-1  ^-1 

F.f- 

-2np fl  (i+?-Tl-^)r|jJdI]d? 
-1  -1 

Fa,  1 

-  -2ttpx  f  f  1  (l+§+T]+?ll)r| jjdT)d§ 
-1-1 

F*.  1 

-  -2tt  Pi  f1  f  (l-?+Tl-P,T1)rjj|dTld? 
-1  *-l 

4.  STRESSES  AT  ELEMENT  CENTROID 
From  earlier  developments  of  strain-displacement  relations  it  is 


Let  p,  be  the  density  of  body  of  element, 
to  dead  load  is  given  by 

n  a 


or 


where 


Fm  *  2rr  a  *^|j|rdT)d5 
,  JTP|j|rdlld? 

\  o  ♦, 


F  ■  2n 


0  ♦.  0 


0  ♦, 


*9  0 


P-  Co  nJ' 


possible  to  write 


150 


v  f~  ~  ~  T 

~  “  |6z  Or  r  Or  OzJ 

or  (Ara  +  Bft5  +  Cr  aT))va 

^  (A**  +  Bt  +  Cx„i|)ua 

v  -  jf^n 

r  #  4C^'(l  +  Bn  5  +  C„  D  +  D„  §T»ua 

Arn  +  Bra?  +  CP  aT])uB  +  (AI#+  Bta  +  Ct#Tl)va 

in  which  all  notations  are  defined  earlier. 

If  stresses  are  sought  at  the  centroid  of  an  element  we  have 
5  *  T)  ■  0,  Therefore, 


The  stresses  at  the  centroid  are  then  obtained  from  the  standard 
stress-strain  relations. 
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APPENDIX  B 
COMPUTER  PROCRAMS 


There  are  six  independent  computer  programs  in  the  present  report: 

(1)  SP1  -  This  program  performs  static  elastoplastic  analysis  for 
a  fiber- reinforced  thin  shell. 

(2)  SVP1  -  This  program  performs  static  viscoelastoplastic  analysis 
for  a  fiber- reinforced  thin  shell. 

(3)  DVP1  -  This  program  performs  dynamic  viscoelastoplastic  analysis 
for  a  fiber-reinforced  thin  shell.  Dynamic  elastoplastic  analy¬ 
sis  is  obtained  as  a  special  case. 

(4)  SF2  -  This  program  performs  static  elastoplastic. analysis  for 
a  fiber- reinforced  thick  shell. 

(5)  SVP2  -  This  program  performs  static  viscoelastoplastic  analysis 
for  a  fiber- reinforced  thick  shell. 

(6)  DVP2  -  This  program  performs  dynamic  viscoelastoplastic  analysis 
for  a  fiber- reinforced  thick  shell.  Dynamic  elastoplastic  analy¬ 
sis  is  obtained  as  a  special  case.  This  program  is  also  capable 
of  analyzing  a  body  composed  of  several  monolithic  and  composite 
materials . 

It  should  be  noted  that  all  of  these  6  programs  can  handle  isotropic  ma¬ 
terials  and  linear  elastic  behavior  as  special  cases. 

A  summary  of  the  basic  theories  used  and  various  features  available 


in  the  programs  is  given  below: 
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1.  Thin  Shell 

a.  Novozhilov’s  linear  thin  shell  theory  is  used. 

b.  Curved  finite  element  is  used.  Interpolation  functions  are 
based  on  linear  variations  of  meridional  and  tangential  dis¬ 
placement  and  cubic  variation  of  transverse  displacement. 

c.  Modified  Hill's  anisotropic  yield  parameters  are  derived  and 
specialized  for  layers  of  angle  or  cross  plys. 

d.  Internal  or  hidden  variables  are  used  for  viscous  effects. 

e.  Solution  of  governing  equations  is  based  on  incremental  approach 
with  increments  in  loading  for  time- Independent  problems  and 
increments  in  time  for  time  dependent  problems. 

2.  Thick  Shell 

a.  Three  dimensional  theory  for  strain-displacement  is  used  and 
specialized  for  axisymmetric  deformations. 

b.  Plane  strain  isoparametric  finite  element  is  used  with  linear 
variation  of  radial  and  axial  displacements. 

c.  Treatment  of  inelastic  and  viscous  properties  and  solution 
procedures  are  the  same  as  in  a  thin  shell. 

Descriptions  of  each  program  including  subroutine  organization  chart, 
subroutine  descriptions,  flow  chart,  and  data  input  format  are  given  in 
the  following  subappendices. 


APPENDIX  B.1.1 


8P1  -  SUBROUTINE  ORGANIZATION  CHART 


— Ibcond  1 

— MD 

HgAftoL  | - j  Speedy  I 

— (GMETRy  i 
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SUBROUTINE 

NAME 

ASSEMB 

BCOND 

DISPL 

DLOAD 

ELIMIT 

ELSTIF 

FACTOR 

SOLTN 

GMETRY 

MATRL 

PLASTC 

PUSH 

RITEXK 

SPEEDY 

STRESS 

ZERO 

ZRDISP 


APPENDIX  B.1.2 

SP1  -  DESCRIPTIONS  OF  SUBROUTINES 


DESCRIPTION 


Assembles  stiffness  matrix  and  force  vector  in  a  single  sub¬ 
scripted  array  with  boundary  conditions  imposed. 

Reads  boundary  conditions  and  initializes  the  row-column 
indices  for  assembling. 

Calls  FACTOR  and  SOLTN,  and  writes  nodal  forces  and  displace¬ 
ments  . 

Computes  the  equivalent  nodal  forces.  '  ■ 


Computes  and  writes  the  elastic  limit  displacements,  elastic 
limit  stresses,  and  elastic  limit  load. 

Forms  elastic  stiffness  matrix. 

Solves  given  simultaneous  equations  in  a  single  subscripted 
array. 

Calculates,  for  each  element,  the  geometrical  quantities 
necessary  for  Simpson's  integration  along  the  line  element. 

Reads  all  material  properties  and  constructs  transformation 
matrices  and  local  and  global  elasticity  matrices. 

Performs  incremental  clastoplastlc  analysis. 

Performs  Simpson's  integration  with  given  stress-strain  rela¬ 
tions  to  form  the  element  stiffness  matrix. 

Writes  non- zero  elements  of  the  lower  half  of  the  global 
stiffness  matrix  with  row  number. 

Performs  matrix  multiplication. 

Computes  strains  and  stresses. 

Zeroes  out  the  given  matrix. 

Transforms  displacements  into  Z-R  coordinate  system  and 


writes. 
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APPENDIX  B.1.3 
SP1  -  FLOW  CHART 


Compute  global  and  local  stresses,  and  write 


NINCR  *  0 


Compute : 

(1)  elastic  limit  stress  (§) 

(2)  elastic  limit  displace¬ 
ment  (0O) 

(3)  incremental  load  - 
parameter  DLOFC 

(4)  incremental  load  (dF) 

Let  FTOT  =  XL 


XL  =  ^'T  /  o  max 
DLOFC  =  XL/DELL 


,  <f5\  =  a /DELL  i  ®  =  ,  a  =  <j> 


DO  NInl,  NINCR 


FTOT  =  FTOT  +  DLOFC 
where  FTOT  is  the  load  parameter 


KOUNT  r-  0 


KOUNT  •  KOUNT  +  1 
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APPENDIX  B.1.4 
SP1  -  DATA  INPUT  FORMAT 

Card  Is  FORMAT  (20A4) 

TITLE  (I)  Title  of  the  problem 

Card  2:  FORMAT  (1015,  3F10.4) 

(1)  NELEMS  -  No.  of  elements 

(2)  NNODES  -  No.  of  nodes 

(3)  NET  -  No.  of  stations  for  Simpson's  integration 

Suggested  values  of  NET 
0b  -  01  NET 

0  -  3°  15 

3°  -  5°  19 

5°  -  9°  23 

9°  -  14°  29 

The  program  assumes  NET  -  15,  if  NET  a  0 

(4)  MECH  -  Signal  for  distributed  load 

0,  if  a  force  vector  for  distributed  loads  is  not  wanted 
1,  if  a  force  vector  for  distributed  loads  Is  wanted 

(5)  NDITO  -  Signal  for  uniform  or  Irregular  distributed  load 

0,  uniform  pressure 
1,  pressure  varies  meridionally 

(6)  NINCR  -  No.  of  load  increment  desired 

Set  NINCR  =  DELL  if  two  times  the  elastic  limit  load 


is  to  be  applied. 
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Suggested  values  of  NINCR  and  DELL 
NINCR  =  40  ~  100 

if  NCY  =  0 

DELL  r,  40  -  100 
NINCR  -  20  ~  50 

if  NCY  *  3 

DELL  -  20  -  50 

NINCR  -  0  for  elastic  analysis 

(7)  NLA  -  No.  of  layers  for  angle  or  cross  plys 

(8)  NBC  -  No.  of  nodes  with  boundary  condition(s) 

(9)  NCY  -  No.  of  iterations  within  a  load  increment 

(10)  NCON  -  No.  of  node(s)  with  concentrated  load (a) 

(11)  DELL  -  Fraction  of  elastic  limit  load  to  be  applied  for  each 

load  increment. 

(12)  PER  -  Percent  error  allowed  for  convergence  in  elastoplastlc 

analysis.  Suggested  value:  5  —  10 

(13)  PCTARC-  Fraction  of  arc  length  ignored  at  a  pole. 

Suggested  value:  .01 
Cards  3:  FORMAT  (3X,  2E12.6) 

(1)  Z(N)  -  Z  -  coordinate  value 

(2)  R(N)  -  R  -  coordinate  value 

Note:  Provide  one  card  for  each  node  in  the  order  of  node  number. 
The  node  number  may  be  recorded  in  3X  spaces. 

Cards  4:  FORMAT  (213,  4E12.6) 

(1)  N0DE1  -  first  node  of  the  element 

(2)  N0DE2  -  second  node  of  the  element 

(3)  PHI1  -  ^ 

(4)  PHI2  -  0, 

(5)  HI  -  thickness  at  N0DE1 
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(6)  H2  -  thickness  at  NODE2 

Note:  Repeat  NELEMS  times  in  the  order  of  element  number. 
The  sign  convention  for  0  is  shown  below 


»>  r 


(2)  NDIR(l) 


(3)  NDIR(2)  = 


(4)  NDIR(3)  = 


(5)  NDIR(4) 


Card(s)  5:  FORMAT  (515) 

(1)  NOD  -  node  number  with  boundary  condltion(s) 

'1,  if  meridional  displacement  is  not  allowed 
JO,  if  frbe 

1,  if  circumferential  displacement  is  not  allowed 
.0,  if  free 

1,  if  normal  displacement  is  not  allowed 
.0,  if  free 

'l,  if  the  rotation  of  the  normal  is  not  allowed 
.0,  if  free 

Note:  Repeat  NBC  times. 

Card  6:  FORMAT  (15) 

'0,  if  fiber  angles  are  same  for  every  element 

IANG  = 

.1,  if  not 

Card  7:  FORMAT  (8F10.0) 

SIGS(I)  -  Yidld  stress  in  the  direction  of  the  normal  to  fiber  for 


each  layer. 


Card  8:  FORMAT  (8F10.0) 

SIGT(I)  -  Yield  stress  in  fiber  direction  for  each  layer. 

Card  9:  FORMAT  (8F10.0) 

SIC4(I)  -  Yield  stress  in  45°  from  fiber  direction  for  each  layer. 
Card  10:  FORMAT  (8F10.0) 

TAUT(I)  -  Yield  stress  in  shear  for  each  layer 
Card  11:  FORMAT  (8F10.0) 

EPX(I)  -  Elastic  modulus  normal  to  fiber  for  each  layer 
Card  12:  FORMAT  (8F10.0) 

EFY(I)  -  Elastic  modulus  in  fiber  direction  for  each  layer 
Note:  EFX(I)  =  EPY(I)  for  isotropic  case 
Card  13:  FORMAT  (8F10.0) 

GP(I)  -  Elastic  shear  modulus  for  each  layer 

Card  14:  FORMAT  (8F10.0) 

XNU(I)  -  Poisson's  ratio  (normal  to  fiber)  for  each  layer 
Card  15:  FORMAT  (8F10.0) 

YNU(I)  -  Poisson's  ratio  (in  fiber  direction)  for  each  layer 
Note:  XNU(I)  =  YNU(I)  for  isotropic  case 
Card  16:  FORMAT  (8F10.0) 

EPX(I)  -  Plastic  modulus  normal  to  fiber  for  each  layer 
Card  17:  FORMAT  (8F10.0) 

EPY(I)  -  Plastic  modulus  in  fiber  direction  for  each  layer 
Card  18:  FORMAT  (8F10.0) 

EP4(I)  -  Flastic  modulus  in  45°  from  fiber  direction  for  each 
layer 


Card  19:  FORMAT  (8F10.0) 

GP(I)  -  Plastic  shear  modulus  for  each  layer 
Card(s)  20:  FORMAT  (8F10.0) 

ALPHA(I)  -  fiber  angles  for  each  layer  measured  from  horlsontal 
line. 

Note:  Repeat  Card(s)  20  NELEMS  times  in  the  order  of 
element  number  if  LANG  4  0. 

Card (s)  21:  FORMAT  (6E10.0) 

(1)  PPl(l)  -  distributed  pressure  in  direction  1  at  node  1 

(2)  PP2(1)  -  distributed  pressure  in  direction  2  at  node  1 

(3)  PP3(1)  -  distributed  pressure  in  direction  3  at  node  1 

(4)  PP1(2)  -  distributed  pressure  in  direction  1  at  node  2 

(5)  PP2(2)  -  distributed  pressure  in  direction  2  at  node  2 

(6)  PP3(3)  «  distributed  pressure  in  direction  3  at  node  2 

Note:  Repeat  NELEMS  times  in  the  order  of  element 
number  if  NDITO  =  1 
Card(s)  22:  FORMAT  (15,  4E15.6) 

(1)  NOD  -  node  number 

(2)  CL(1)  -  concentrated  force  in  direction  1 

(3)  CL(2)  -  concentrated  force  in  direction  2 

(4)  CL(3)  -  concentrated  force  in  direction  3 

(5)  CL(4)  -  concentrated  force  in  direction  4 

Note:  Repeat  NCOK  times. 


Omit  if  NCON  =  0. 


mm® 


164 


APPENDIX  B.2.2 

SVP1  -  DESCRIPTIONS  OF  SUBROUTINES 

SUBROUTINE 

NAME  DESCRIPTION 

(r)  (  r)  ( r  ) 

ABG  Calculates  viscous  constants  A,  B,  C. 

ASBLYV  Assembles  in  global  form. 

BCONDY  Imposes  boundary  condition(s). 

CMAT  Forms  the  viscoaity  matrix. 

DLOAD  Computes  the  equivalent  nodal  forces. 

ELSTIF  Reads  the  material  properties  and  forms  the  elastic  stiff¬ 

ness  matrix. 

FORV  Performs  integration  for  viscous  force  vector. 

GMETRY  Calculates,  for  each  element,  the  geometrical  quantities 

necessary  for  Simpson's  integration  along  the  line  element. 

INIT  Reads  the  viscous  property  and  forms  submatrices  for  vis¬ 

cosity  matrix. 

OUTDIS  Writes  the  displacements. 

OUTSTR  Writes  the  stresses  and  the  equivalent  stresses. 

PLASTI  Performs  the  Simpson's  integration  to  form  stiffness  or 

viscosity  matrix. 

PLASV1  Develops  plastic  stiffness  matrix. 

PLASV2  Develops  plastic  stiffness  matrix. 

SPEEDY  Performs  matrix  multiplication. 

STAN  Calculates  viscous  and  elastic  stresses,  checks  yielding. 

STRAW  Calculates  Incremental  viscous  and  plastic  stresses,  in¬ 

cremental  equivalent  stresses,  and  per  cent  error. 

VISCP  Performs  step  by  step  integration. 
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APPENDIX  B.2.2  (Cont.) 

SUBROUTINE 

NAME  DESCRIPTION 

XKINVS  Inverts  a  given  symmetric  matrix. 

ZERO  Zeroes  out  a  given  matrix. 

ZRDISP  Transforms  the  displacements  into  the  Z-R  coordinate 

system  and  writes. 
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APPENDIX  B.2.4 

SVP1  -  DATA  INPUT  FORMAT 

Card  1:  FORMAT  (20AA) 

TITLE (I)  Title  of  the  Problem 
Card  2:  FORMAT  (1015,  3F10.4) 

(1)  NELEMS  -  No.  of  elements 

(2)  NNODES  -  No.  of  nodes 

(3)  NET  -  No.  of  stations  for  Simpson’s  integration 

Suggested  vr .ue  of  NET: 

tp2  -  NET 

0-3°  15 

3°  -  5°  19 

5°  ~  9°  23 

9°  -  15°  29 

The  program  assumes  NET  “15,  if  NET  •  0. 

(4)  MECH  -  Signal  for  distributed  load 

0,  if  a  force  vector  for  distributed  load  is  not 
wanted . 

1,  if  a  force  vector  for  distributed  load  is  wanted. 

(5)  NDITO  -  Signal  for  uniform  of  irregular  distributed  load 

0,  uniform  pressure 
1,  pressure  varies  meridionally 


(6)  NDELT  -  No.  of  time  steps  desired. 

(7)  NLA  -  No.  of  layers  for  angle  or  cross  plys. 

(8)  NBC  -  No.  of  boundary  condition(s). 
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(9)  NCY  -  No.  of  iterations  within  a  time  step. 

(10)  NCCN  -  No.  of  node(s)  with  concentrated  load(a). 

(11)  PER  -  Percent  error  allowed  for  convergence  in  plastic 

analysis. 

Suggested  value:  PER  *  5  ~  10 

(12)  PCTARC  -  Fraction  of  arc  length  ignored  at  a  pole. 

Suggested  value:  0.01 

The  program  assumes  PCTARC  *  0.01,  if  PCTARC  ■  0. 


Card  3: 

FORMAT  (415) 

(1) 

IS(1) 

-  1, 

if  elastic 

analysis  is  wanted. 

0, 

if  plastic 

analysis  is  wanted 

(2) 

IS(2) 

-  1, 

if  nonviscous  analysis  is  wanted 

-  0, 

if  viscous 

analysis  is  wanted 

(3) 

IS(3) 

-  Print  signal 

til 

Print  displacements  and  stresses  at  every  IS(3) — 
step. 

Tae  program  assumes  IS(3)  *  5,  if  IS(3)  "  0. 

(4)  IS(4)  -  Transformation  sighal. 

Transform  the  displacements  into  the  Z-R  coordi- 

til 

nate  system  and  print  for  every  IS(4) —  Btep, 

The  program  assumes  IS(4)  ■  NDELT,  if  IS(4)  "  0. 

Card  4:  FORMAT  (3X,  E12.6) 

DELT  -  Size  of  the  time  step  in  seconds  (At). 

Cards  5:  FORMAT  (3X,  2E12.6) 


(1)  Z(N)  -  Z  *  coordinate  value 

(2)  R(N)  -  R  -  coordinate  value 


Mote:  Provide  one  card  for  each  node  in  the  order  of  node 
number.  The  node  number  may  be  recorded  in  3X  spaces. 
Cards  6:  FORMAT  (213,  4E12.6) 

(1)  N0DE1  -  First  node  of  the  element 

(2)  N0DE2  -  Second  node  of  the  element 

(3)  PHI1  -  «PX 

(4)  PHI 2  -  <Pa 

(5)  HI  -  Thickness  at  N0DE1 

(6)  H2  -  Thickness  at  N0DE2 

Mote:  Repeat  NELEMS  times  in  the  order  of  element  number 
The  sign  convention  for  <P  is  shown  below. 


Card(s)  7:  FORMAT  (6E10.0) 

(1)  PP1(1)  -  Distributed  pressure  in  direction  1  at  node  1. 

(2)  PP2(1)  -  Distributed  pressure  in  direction  2  at  node  1. 

(3)  PP3(1)  -  Distributed  pressure  in  direction  3  at  node  1. 

(4)  PP1(2)  -  Distributed  pressure  in  direction  1  at  node  2. 

(5)  PP2(2)  -  Distributed  pressure  in  direction  2  at  node  2. 

(6)  PP3(3)  -  Distributed  pressure  in  direction  3  at  node  2. 
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Note:  Repeat  NELEMS  times  in  the  order  of  element  number  if 
NDITO  -  1. 

Card(s)  8i  FORMAT  (15,  4E15.6) 

(1)  NOD  -  Node  number 

(2)  CL(1)  -  Concentrated  force  in  direction  1. 

(3)  CL(2)  -  Concentrated  force  in  direction  2. 

(4)  CL(3)  -  Concentrated  force  in  direction  3. 

* 

(5)  CL(4)  -  Concentrated  force  in  direction  4, 

Note:  Repeat  NCON  times. 

Omit  if  NCON  «  0. 

Card(s)  9:  FORMAT  (215) 

(1)  NOD  -  Node  number  with  boundary  condition(s), 

(2)  NDEGF  -  Thi  coordinate  number  of  whicn  freedom  is  re¬ 

strained. 

Note:  Repeat  NBC  times. 

Card  10:  FORMAT  (8F10.0) 

SIGS(I)  -  Yield  stress  in  the  direction  of  the  normal  to  fiber 
for  each  layer. 

Card  11:  FORMAT  (8F10.0) 

S"!GT(I)  -  Yield  stress  .  i  fiber  direct*  for  each  layer. 

Card  12:  FORMAT  (8F10.0) 

SIG4(I)  -  Yield  stress  in  45°  from  fiber  direction  for  each  laye*. 
Card  13:  FORMAT  (8F10.0) 

TAUT(l)  -  Yield  stress  in  shear  for  each  layer. 

Card  14:  FORMAT  (8F10.U) 

EPX(l)  -  Elastic  modulus  normal  to  fiber  for  each  layer. 
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Card  15:  FORMAT(8FlO.O) 

EPY(I)  -  Elastic  modulus  in  fiber  direction  for  each  layer. 

Ntte:  EPX(I)  •  EPY(I)  for  isotropic  case. 

Card  16:  FORMAT  (8F10.0) 

GP(I)  -  Elastic  shear  modulus  for  each  layer. 

Card  17:  FORMAT  (8F10.0) 

XNU(I)  -  Poisson's  ratio  (normal  to  fiber)  for  each  layer. 

Card  18:  FORMAT  (SF10.0) 

TNU(I)  -  Poisson's  ratio  (in  fiber  direction)  for  each  layer. 
Note:  XNU(I)  *  YNU(I)  for  isotropic  case. 

Card  19:  FORMAT  (8F1.0.0) 

EPX(I)  -  Plastic  modulus  normal  to  fiber  for  each  layer. 

Card  20:  FORMAT  (8F10.0) 

EPY(I)  -  Plastic  modulus  in  fiber  direction  for  each  layer. 

Card  21:  FORMAT  (8F10.0) 

EPr(l)  -  Plastic  modulus  in  45°  from  fiber  direction  for  each 
layer. 

Card  22:  FORMAT  (8F10.0) 

GP(I)  -  Plastic  shear  modulus  for  each  layer. 

Card  23:  FORMAT  (8F10.0) 

ALPHA(I)  -  Fiber  angles  for  each  layer  measured  from  horizontal 
line. 

Note:  Repeat  NELEMS  times. 

Card  24:  FORMAT  (8F1C.0) 

(1)  r;.  -  Relaxation  time  in  fiber  direction. 

(2)  TT  -  Relaxation  time  in  the  direction  of  normal  to  fiber. 


173 


(3) 

(4) 
Note: 


DSL  -  EPY(I)  (see  Card  15) 
DST  -  EPX(I)  (aee  Card  14) 
Not  required  if  IS(2)  j1  0. 
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APPENDIX  B.3.1 

DVP1  -  SUBROUTINE  ORGANIZATION  CHAPT 
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APPENDIX  B.3.2 

DVP1  -  DESCRIPTIONS  OF  SUBROUTINES 


SUBROUTINE 

NAME  DESCRIPTION 


ABG 

^r)(r)lr) 

Calculates  viscous  constants  a,  b»  C* 

ASBKTV 

Assembles  In  global  form. 

BCONDY 

Imposes  boundary  condltlon(s). 

CMAT 

Forms  the  viscosity  matrix. 

DLOAD 

Computes  the  equivalent  nodal  forces 

ELSTIF 

Reads  the  material  properties  and  forms  the  elastic  stiff¬ 
ness  matrix. 

FORV 

Performs  integration  for  viscous  force  vector. 

GMETRY 

Calculates,  for  each  element,  the  geometrical  quantities 
necessary  for  Simpson's  integration  along  the  line  ele¬ 
ment. 

INIT 

Reads  the  viscous  property  and  forms  submatrices  for 
viscosity  matrix. 

MASS 

Forms,  for  each  element,  the  consistent  mass  matrix. 

OUTDIS 

Writes  the  displacements. 

OUTSTR 

Writes  the  stresses  and  the  equivalent  stresses. 

PLASTI 

Performs  Simpson's  integration  to  form  stiffness  or 
viscousity  matrix. 

PLASV1 

Develops  plastic  stiffness  matrix. 

PLASV2 

Develops  plastic  stiffness  matrix. 

SPEEDY 

Performs  matrix  multiplication. 

STAN 

Calculates  viscous  and  elastic  stresses,  checks  yielding. 

STRAUP 

Calculates  incremental  viscous  and  plastic  stresses, 
incremental  equivalent  stresses,  and  per  cent  error. 

VISCP 

Performs  step-bystep  integration. 
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APPENDIX  B.3.2  (cent.) 

DVPI  -  DESCRIPTIONS  OF  SUBROUTINES 

SUBROUTINE 

NAME  DESCRIPTION 

XKINVS  Inverts  a  given  symmetric  matrix. 

ZERO  Zeroes  out  a  given  matrix. 

Transforms  the  displacements  into  the  Z-R  coordinate 
system  and  writes. 


ZRDISP 


APPENDIX  B.3.3  (continued) 
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APPENDIX  B.3.4 
DVP1  -  DATA  INPUT  FORMAT 


Card  1:  FORMAT  (20A4) 

TITLE (1)  Title  of  the  problem 

Card  2s  FORMAT  (1015,3F10.4) 

(1)  NELEMS  -  No.  of  elements 

(2)  NNODES  -  No.' of  nodes 

(3)  NET  -  No.  of  stations  for  Simpson's  integration 

Suggested  value  of  NET 


eft,-  <Px 

NET 

0  ~  3° 

15 

3°  ~  5° 

19 

5°  -9° 

23 

9°  ~  15° 

29 

The  program  assumes  NET  "15,  if  NET  ■  0. 

(4)  MECH  -  Signal  for  distributed  load 

0,  if  a  force  vector  for  distributed  load  is  not 
wanted . 

1,  if  a  force  vector  for  distributed  load  is  wanted. 

(5)  NDITO  -  Signal  for  uniform  of  irregular  distributed  load 

0,  uniform  pressure 
1,  pressure  varies  meridionally 

(6)  NDELT  -  No.  of  Htne  steps  desired 

(7)  NLA  -  No.  of  layers  for  angle  or  cross  plys 

(8)  NBC  -  No.  of  boundary  condltlon(s) 


1 


Transform  Che  displacements  into  Che  r-R  coordinate 
system  and  print  for  every  IS(4)^  step. 

The  program  assumes  IS(4)  •  NDELT,  if  IS(4)  ■  0. 
Card  4:  F0RMAT(3X,2E12.6) 

(1)  RHO  -  Average  density  of  the  material  to  compute  the  mass 

matrix  In  Ibs/irt* 

(2)  CELT  -  Sire  of  the  time  step  in  seconds  (At) 

Note:  At  *  10"8  sec.  for  metals. 
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| 

< 

(9) 

NCY 

-  No.  of  iterations  within  a  time  step  1 

"i 

i 

(10) 

NCON 

-  No.  of  node(8)  with  concentrated  load(s)  { 

(ID 

PER 

-  Percent  error  allowed  for  convergence  in  plastic 

analysis  j 

< 

\ 

Suggested  value  PER  *  5  ~  10  j 

i 

* 

3 

J 

(12) 

PCTARC 

-  Fraction  of  arc  length  ignored  at  a  pole 

j 

i 

Suggested  value:  0.01  j 

* 

The  program  assumes  PCTARC  “  0.01,  if  PCTARC  “  0. 

l 

A 

Card  3: 

FORMAT  (415)  j 

i 

/} 

‘{ 

\ 

y 

i 

5 

1 

(1) 

IS(1) 

* 

-  1,  if  elastic  analysis  is  wanted  j 

| 

0,  if  plastic  analysis  is  wanted  ; 

* 

i 

i 

(2) 

IS(2) 

-  1,  if  nonvlscout  analysis  is  wanted 

& 

■■t 

j 

0,  if  viscous  analysis  is  wanted 

1 

(3) 

IS(3) 

-  Print  signal 

Print  displacements  end  stresses  at  every  IS  (3)— 

I 

* 

i 

5 

time  step. 

; 

The  program  assumes  IS(3)  ■  5,  if  IS(3)  ■  0 

(4) 

IS(4) 

••  Transformation  signal 

Cards  5:  FORMAT  (3X.2F12.6) 

(1)  Z(H)  -  Z  -  coordinate  value 

(2)  R(N)  -  R  -  coordinate  value 

Note:  Provide  one  card  for  each  node  in  the  order  of  node  number. 
The  node  number  may  be  recorded  in  3X  spaces. 

Cards  6:  FORMAT  (213,  4E12.6) 

(1)  N0DE1  -  First  node  of  the  element 

(2)  NODE2  -  Second  node  of  the  element 

(3)  PHI1  - 

(4)  PHI2  -  qx, 

(5)  HI  -  Thickness  at  N0DE1 

(6)  H2  -  Thickness  at  NODE2 

Note:  Repeat  NELEMS  times  in  the  order  of  element  number. 
The  sign  convention  for  <P  is  shown  below. 


Card(s)  7:  FORMAT  (6E10.0) 

(1)  PP1(1)  -  Distributed  pressure  in  direction  1  at  node  1 

(2)  PP2(1)  -  Distributed  pressure  in  direction  2  at  node  1 

(3)  PP3(1)  -  Distributed  pressure  in  direction  3  at  node  1 

(4)  PP1(2)  -  Distributed  pressure  in  direction  1  at  node  2 
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(5)  PP2(2)  -  Distributed  pressure  In  direction  2  at  node  2 

(6)  PP3(3)  -  Distributed  pressure  In  direction  3  at  node  2 

Note:  Repeat  NELEMS  times  In  the  order  of  element 
number  if  NDITO  ■  1. 


Card(s) 

8:  FORMAT  (15,  4E15.6) 

(1) 

NOD 

-  Node  number 

(2) 

CL(1) 

-  Concentrated 

force 

in  direction 

(3) 

CL(2) 

-  Concentrated 

force 

in  direction 

(4) 

CL(3) 

-  Concentrated 

force 

in  direction 

(5) 

CL(4) 

-  Concentrated 

force 

in  direction 

Note:  Repeat  NOON  times. 

Omit  if  NOON  -  0. 

Card(s)  9:  FORMAT  (215) 

(1)  NOD  '  Node  number  with  boundary  condltion(s) 

(2)  NDEGF  -  The  coordinate  number  of  which  freedom  Is  restrained. 

Note:  Repeat  NBC  times. 

Card  10:  FORMAT  (8F10.0) 

S1GS(I)  -  Yield  stress  in  the  direction  of  the  norsial  to 
fiber  for  each  layer. 

Card  11:  FORMAT  (8F10.0) 

SIGT(I)  -  Yield  stress  in  fiber  direction  for  each  layer. 

Card  12:  FORMAT  (8F10.0) 

SIGA(I)  -  Yield  stress  in  45°  from  fiber  direction  for  each 
layer 

Card  13:  FORMAr  (8F10.0) 

TAUT(I)  *  Yield  stress  in  shear  for  each  layer 


|  i 
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Card  14:  FORMAf  (8F10.0) 

BPX(I)  *  Elastic  modulus  normal  to  fiber  for  each  layer 
Card  15:  FORMAT  (8F10.0) 

EPY(I)  -  Elastic  modulus  in  fiber  direction  for  each  layer 
Note:  FPX(I)  ■  EPY(I)  for  isotropic  case 
Card  16:  FORMAT  (8F10.0) 

GP(I)  -  Elastic  shear  modulus  for  each  layer 
Card  17:  FORMAT  (8F10.0) 

XNU(I)  -  Poisson's  ratio  (normal  to  fiber)  for  each  layer 
Card  18:  FORMAT  (8F10.0) 

YHU(I)  -  Poisson's  ratio  (in  fiber  direction)  for  each  layer 
Mote:  XNU(I)  *  YNU(I)  for  isotropic  case 
Card  19:  FORMAT  (8F10.0) 

EPX(l)  -  Plastic  modulus  normal  to  fiber  for  each  layer 
Card  20:  FORMAT  (8F10.0) 

EPY(I)  -  Plastic  modulus  in  fiber  direction  for  each  layer 
Card  21:  FORMAT  (8F10.0) 

EP4(I)  -  Plastic  modulus  in  45°  from  fiber  direction  for  each 
layer 

Card  22:  FORMAT  (EriO.Q) 

GP(I)  •  Plastic  shear  modulus  for  each  layer 
Card  23:  FORMAT  (8F10.0) 

ALPHA(I)  -  Fiber  angles  for  each  Layer  measured  from  horlzonta 
line 

Mote:  Repeat  MEhEMS  times. 


Card  24:  FORMAT  (8F10.0) 


(1)  TL  -  Relaxation  time  In  fiber  direction 

(2)  TT  -  Relaxation  time  in  the  direction  of  normal  to  fiber 

(3)  DSL  -  EPY(I)  (see  Card  15) 


APPENDIX  B.4.2 


SUBROUTINE 

NAME 

ASBLYI 

F 

GAUSS 

INPUT 

PUSH 

PUSI2 

SETUP 

SPEEDY 

STIFIP 

STIFF 1 

STRANP 

STRESS 

TRANS 

XKINVS 

XTIFIS 


SP2  -  DESCRIPTION  OF  SUBROUTINES 


DESCRIPTION 


Assembles  stiffness  matrix. 

Evaluates  functional  values  used  in  Gaussian  quadrature 
integration. 

Performs  Gaussian  quadrature  integration. 

Reads  and  writes  input  data. 

Develops  instantaneous  stiffness  matrix  for  initial 
yielding. 

Develops  instantaneous  stiffness  matrix  for  subsequent 
yielding. 

Sets  up  weight-function  values  used  in  quadrature  inte* 
gration. 

Performs  matrix  multiplication. 

Forms  element  stiffness  matrix  for  the  given  elasticity 
matrix. 

Calculates  coordinate-dependent  values  used  in  quadrature 
Integration. 

Calculates  incremental  stresses,  incremental  equivalent 
stresses,  and  ner  cent  error  to  test  for  convergence. 

Computes  displacements  and  stresses  under  given  load, 
and  then  scales  the  quantities  to  the  elastic  limit. 

Develops  transformation  matrices  for  both  stress  and  strain. 

Performs  matrix  inversion. 

Calculates  local  elasticity  matrix,  transforms  this  to 
global  coordinates,  and  then  forms  global  stiffness  matrix. 
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APPENDIX  B.4.3 
SP2  -  FLOW  CHART 


188 


_  WATT\TfP  1  1 


dJT  -  B  d® 


da  -  (D  +  D)dE 


da8  «  f (do ) 


Check  Convergence 


- <T  KOUNT  ■  :  NCYJ> - Jdax  -  da3 


Update  previous  quantities,  1,  e. 

0j  -  ®,  +  d®, 

a i  ■  a i  -j  +  dOj 

*  a(  -i  +  do 


<; 


STOP 


i}  ' 


-  iV-' 


Card  1: 


Card  2: 


APPENDIX  B.4.4 
SP2  -  DATA  INPUT  FORMAT 


FORMAT  (20A4) 

TITLE (I)  Title  of  the  problem 
FORMAT  (215,  2F10.0) 


(1)  NCY  -  Number  of  iterations  allowed  for  convergence  within 

a  load  increment, 

(2)  NINCR  -  Number  of  load  increments, 

(3)  PER  -  Allowable  per  cent  error  for  convergence 

(4)  DELL  -  Fraction  of  elastic  limit  load  to  be  applied  for 

each  load  increment. 

Cards  3:  FORMAT  (2F10.0) 

(1)  z  -  z  coordinate  value 


(2)  R 


-  R  coordinate  vlaue 


(Note:  Provide  one  card  for  each  node. 

Cards  4:  FORMAT  (415) 

I-J-K-L  designations  of  nodes  for  a  given  element,  in  counter¬ 
clockwise  order. 

Note:  Provide  one  card  for  each  element. 


Sample: 


1-4 
J  -  5 
K  -  2 
L  -  1 


Card  5:  FORMAT  (15) 

NBC  -  Number  of  boundary  conditions 

Cards  6:  FORMAT  '215) 


(1)  NOD  -  Node  number  of  boundary  condition 

(2)  IDG  -  Degree  of  freedom  at  the  node  which  is  restrained 
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IDG  *  l  if  z  displacement  restrained 
IDG  '•*  2  if  R  displacement  restrained 
Card  7:  FORMAT  (15) 

ISAME  -  ISAME  *  1  if  all  elements  have  same  card(s)  8. 

Card(s)  8:  FORMAT  (8F10.0) 

(1)  ET  -  Modulus  of  elasticity  transverse  to  fiber 

(2)  EL  -  Modulus  of  elasticity  along  fiber 

(3)  GNUTL  -  Poisson's  ratio  between  fiber  and  direction  trans¬ 

verse  to  fiber 

(4)  GNUTT  -  Poisson's  ratio  between  directions  normal  to  fiber 

(5)  GTL  -  Shear  modulus  between  fiber  and  direction  trans¬ 

verse  to  fiber 

(6)  GT  -  Shear  modulus  between  direction*  normal  to  fiber 

(7)  ALPH  -  Angle  in  degrees  from  horizontal  to  fiber  direc¬ 

tion  within  plane  of  the  structure 

(8)  PHI  -  Angle  in  degrees  from  vertical  to  plane  of  the 

structure 

Note:  Provide  one  card  for  each  element  if  ISAME  not 
equal  to  1, 

Card  9:  FORMAT  (15) 

ISAME  -  ISAME  *  1  if  each  element  has  same  card(s)  (10)  and  card(s) 

11. 

Card(s)  10:  FORMAT  (8F10.0) 

(1)  SOT  -  Yield  stress  in  tension  transverse  to  fiber 

(2)  SOL  -  Yield  stress  in  tension  in  fiber  direction 

(3)  TOT  -  Shear  yield  stress  in  tension  transverse  to  fiber 

(4)  TOL  -  Shear  yield  stress  in  tension  in  fiber  direction 

(5)  EPT  -  Plastic  modulus  normal  to  fiber 

(6)  EPL  -  Plastic  modulus  along  fiber 


(7)  GPT  -  Plastic  shear  modulus  between  directions  normal 

to  fiber 

(8)  GPL  -  Plastic  shear  modulus  between  fiber  and  direction 

normal  to  fiber 

Card(s)  11:  FORMAT  (8F10.0) 

Same  as  cards  8  except  for  compression 

Note:  Provide  one  card(s)  10  and  one  card(s)  11  for  each 
element  if  ISAME  not  equal  1. 

Card  12:  FORMAT  (15) 

NFORCE  -  Number  of  externally  applied  concentrated  fOrceB 
Card(s)  13:  FORMAT  (215,  F10.0) 


(1)  NCD 


Node  number  of  applied  concentrated  force 


(2)  IDG  -  Degree  of  freedom  at  which  force  is  applied 

IDG  “  1  if  force  applied  in  z  direction 
IDG  ■  2  if  force  applied  in  R  direction 

(3)  FOR  -  Magnitude  of  applied  force 

Noto:  Provide  one  card  for  each  applied  force  if  NFORCE 
=  0  Card  13  is  not  needed 

Card  14:  FORMAT  (15,  F10.0) 

NPRES  -  Number  of  elements  with  applied  pressure 

PRESS  -  Applied  pressure 


Card(s)  15: 


FORMAT  (315) 


(1)  NELP  -  element  pressure  applied  to 

(2)  MODI  -  One  node  of  the  element  having  pressure 

(3)  N0D2  -  Second  node  of  the  element  having  pressure 

Note:  Provide  one  card  for  each  element  having  pressure.  If 
NPRESS  ™  0  card(8)  15  not  needed. 


SAWPLT 


PRESS  «  5.0 
NELP  -  5 
NODI  -  3 
NOD2  -  4 


►.v  /  < 
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APPENDIX  B.5.1 

SVP2  -  SUBROUTINE  ORGAINTZATION  CHART 


-flNPUT  ~*1* 
*  |  SETUP  I 


TRANS 


XTIFIS 


-[  STIFF1 1- 
•I  SPEEDY  I 

-f  stifip] 


(GAUSS 


I INITI 


ASBLYI] 


SPEEDY  I 


STIFIP  | 


I MAIN  I - 


|  XKINVS ] 


I  ASBLYI I 


!  SPEEDY  I 


•I  StAn  I- 


STRAIN] 


SPEEDYl 


SPEEDY 


- [isrer] 


iSFEEDY! 


IPLASI 


- VtokSt  I 


1  STRANPl* 


|  ASBLYI! 


[SPEEDY] 


[SPEEDYj 


■fsPEEDYl 
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SVP2  -  DESCRIPTION  OF  SUBROUTINES 

SUBROUTINE 

NAME _ DESCRIPTION _ 

ABG  Calculates  viscoelastic  constants  A>  B,  C 

ASBLYI  Assembles  stiffness  matrix  and  damping  matrix 

BCFOR  Applies  boundary  conditions  to  force  vectors 

F  Evaluates  functions  used  in  integration 

GAUSS  Performs  Gaussian  quadrature  integration 

INITI  Develops  the  damping  matrix  together  with  values  needed 

in  viscous  stress  calculations 

INPUT  Reads  and  writes  input  data 

PLASH  Develops  plastic  stiffness  matrix  for  initial  yielding 

PLASI2  Develops  plastic  stiffness  matrix  for  subsequent 

yielding 

SPEEDY  Performs  matrix  multiplication 

SETUP  Initializes  values  used  in  Gaussian  quadrature  integration 

STAN  Calculates  elastic  and  viscous  stresses  and  checks  yield* 

ing.  If  no  yielding  it  updates  stresses,  equivalent  stresses 
etc.  If  yielded  it  initializes  Incremental  equivalent  stress 
STIFF1  Calculates  coordinate  values  used  in  the  integration  scheme 

for  development  of  the  stiffness  and  damping  matrix 
STIFIP  Forms  element  stiffness  and  damping  matrix  from  functions 

evaluated  in  STIFF 1 

STRAIN  Calculates  local  strains  from  global  displacements 

Calculates  increment  stresses,  equivalent  stresses  and  per 
cent  error  between  any  two  iterations  of  a  time  increment 


STRANP 


» *Trr*v 'nwjt "» 


•h  <?  i  >**'■'** ■ ' 
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SUBROUTINE 

NAME _ DESCRIPTION _ _ 

TRANS  Forms  the  stress  and  strain  transformation  matrices 

XKINVS  Performs  matrix  inversion 

XTIFIS  Calculates  and  assembles  global  stiffness  matrix 


APPENDIX  B.5.3 
SVP2  -  FLOW  CHART 


The  flow  chart  for  SVP2  is  the  same  as  that  for  SVP1  in  Appendix 
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APPENDIX  B.5.4 
SVP2  -  DATA  INPUT  FORMAT 

Card  1:  FORMAT  (15,  F10.0) 

(1)  NDELT  -  Number  of  time  increments 

(2)  DELT  -  Size  of  time  increment 

Card  2:  FORMAT  (15,  F10.0) 

(1)  NCY  -  Number  of  iterations  allowed  for  convergence  with¬ 

in  a  time  increment  if  yielding  has  occurred. 

(2)  PER  -  Allowable  per  cent  error  for  convergence 

Cards  3:  Insert  Cards  3  through  Cards  15  from  B.4.4. 

through 

Cards  15 

Card(s)  16: FORMAT  (4F10.0,  15) 

(1)  TT  -  Relaxation  time  in  direction  transverse  to  fiber 

(2)  TL  -  Relaxation  time  in  fiber  direction 

(3)  DT  -  Modulus  of  elasticity  normal  to  fiber 

(4)  DL  -  Modulus  of  elasticity  in  fiber  direction 

(5)  ISAME  -  Let  ISAME  -  1  if  all  elements  have  same  properties 

Note:  Provide  one  card  for  each  element  if  ISAME  not  equal  1. 
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APPENDIX  B.6,1 

DVP2  -  SUBROUTINE  ORGANIZATION  CHART 
The  eubrouClne  organization  chart  for  DVP2  ia  the  same  as  that 
for  SVP2  except  that  MAIN  calls  an  additional  subroutine,  MASS,  which 
develops  the  mass  matrix. 


APPENDIX  B.6.2 

DVP2  -  DESCRIPTION  OF  SUBROUTINES 
The  description  of  subroutines  for  DVP2  is  the  same  as  that  for 
SVP2  except  that  subroutine  MASS  should  be  added  to  the  list. 

MASS  Develops  and  assembles  the  mass  matrix. 


APPENDIX  B.6.3 
DVP2  -  FLOW  CHART 

The  flow  chart  for  DVP2  is  the  same  as  that  for  DVP1  in  Appendix 
B.3.3. 


APPENDIX  B.6.4 
DVP2  -  DATA  INPUT  FORMAT 

Cards  1  -  Insert  Card  1  through  Cards  16  from  SyP2  in  Appendix  B.5.4, 

through 

Card(s)  16 

Card(s)  17:  FORMAT  (F10.0,  15) 

(1)  DEN  -  Weight  density  of  the  material  in  pounds  per 

cubic  inch. 

(2)  ISAME  -  Let  ISAME  ■  1  if  all  elements  have  same  density. 

Note:  Provide  one  card  for  each  element  if  ISAME  not  equal  to  1. 


i 
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